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We present both analytically and numerically the consistent analysis from the Affleck-Dine (AD) 
dynamics to the subsequent semiclassical evolution in both gravity-mediated and gauge-mediated 
models. We obtain analytically the elliptic motions in the AD dynamics as the analogy of the 
well-known Kepler-problem, and by solving the equations of motion in a lattice, we find that the 
0^ ■ semiclassical evolution goes through three distinct stages as a nonequilibrium process of reheating 

' the Universe: pre-thermalisation, bubble collisions and thermalisation. We report that the second 

stage of our case lasts rather long compared to the second stage of the reheating case, and the 
thermalisation process is unique due to the presence of "thermal Q-balls" . 
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INTRODUCTION 



The present baryon asymmetry in the Universe is one of the most mysterious problems in cosmology 
and particle physics, see a review Within the Standard Model (SM), the electroweak baryogen- 
esis was suggested to explain the inequality between baryon and anti-baryon number, and the recent 
_ , developments are shifted into constructing the theory of reheating the Universe [2]. The electroweak 

baryogenesis satisfies the well-known Sakharov's three conditions required for successful baryogenesis Q , 
(— I ' namely baryon number production, charge parity (CP) violation and the process taking place out-of- 

equilibrium; however, the predicted CP violation is too small to explain the present observed baryon 
number. By satisfying the above three conditions, the Affleck-Dine (AD) baryogenesis [3], which was 
proposed in the theoretical framework beyond the SM, namely, the Minimal Super-symmetric Standard 
■ Model (MSSM), is a more successful scenario to tackle this puzzle, since it may solve problems of grav- 

I itino and moduli overproduction and give rise simultaneously to the ordinary matter and dark mater in 

the Universe. The MSSM has many gauge-invariant flat directions along which R parity is preserved. 
The flat directions are lifted by super-symmetry (SUSY) breaking effects arising from nonrenormalisable 
. terms, which give a U(l) violation through A-terms. In the original scenario of the AD baryogenesis, one 

' can parametrise one of the flat directions in terms of a complex scalar field known as an AD field (or 

On , AD condensate which consists of a combination of squarks and/or sleptons fields). The AD field evolves 

' to a large field expectation value during an inflationary epoch in the early Universe. After inflation, the 

orbit of the AD field can be kicked along the phase direction due to the A-terms which generate the 
U(l) charge (baryon/lcpton number), and then the A-terms become negligible, where the AD field rotate 
towards the global minimum of the scalar potential. Hence, the generated global U(l) charge is fixed and 
the orbit of the AD field rotates around the origin of the complex field-space, c. /. the anomaly mediated 
models fs'l. After the AD condensate decays into usual baryons and leptons, AD baryogenesis becomes 
complete. 

The trajectory of the AD field is identical to the planetary orbits in the well-known Kepler-problem 
as we will show later, replacing the Newtonian potential by an isotropic harmonic oscillator potential 
0. This coincidental classical- mechanics reduction was noted for the orbits of a probe brane in the 
Branonium system 0, [!]■ As general relativity predicted that planetary orbits precess by adding the 
relativistic correction to the Newtonian potential, we will see the similar events occur for the orbits of 
AD fields, which are disturbed by quantum and nonrenormalisable effects. 

By including quantum corrections [lo| and/or thermal effects 11 1 in the mass term of the standard 



AD scalar potentials, the AD condensate is classically unstable against spatial perturbations due to the 
jresence of negative pressure [l2| , and fragments to bubble-like objects, eventually evolving into Q-balls 



13 1 . Lee pointed out [ij] that Q-balls may form due to bubble nucleation (first order phase transition) 
15l |. even in the case that the condensate is classically stable against the linear spatial perturbations. 

A Q-ball is a nontopological soliton [l^ whose stability comes from the existence of a continuous global 
or local charge Q (see a review [l3| and references therein). Tsumagari et. al. [H, [l^ showed previously 
the complete stability analysis of Q-balls at zero-temperature in both polynomial potentials and MSSM 
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flat potentials. Laine et. al. [20| investigated the stability of Q-balls in a thermal bath. The stability of 
the thermal SUSY Q-balls is different from the one of the standard "cold" Q-balls, since they suffer from 
evaporation diffusion [2l|, dissociation [lH, and decays into light/massless fermions [Ij]. Therefore, 
most SUSY Q-balls are generally not stable but long-lived, and it may thermalise the Universe by decaying 
into baryons on their surface [24| . which could solve the gravitino and moduli over-production problems 
without fine-tuning. The SUSY Q-balls in gravity-mediated (GRV-M) models are quasi-stable decaying 
into the lightest SUSY particles (LSP dark matter), and the fraction of the baryons from the Q-balls 
may give the present baryon number, which can explain the similarity of the energy density between 
the observed baryons and dark matter [S^, [2^. The SUSY Q-balls in gauge- mediated (GAU-M) models, 
however, can be extremely long-lived so that those Q-balls are a candidate of cold dark matter [2^ and 
may give the present observed baryo-to-photon ratio [l^l • 

The dynamics and formation of Q-balls have been investigated numerically. With different relative 
phases and initial velocities, the authors [26] found a charge transfer from one Q-ball to the other and 
interesting ring formation after the collision. It has been found [131 that similar ring-like solutions 
are responsible for the excited states from the ground state (Q-ball) by introducing extra degrees of 
freedom: spatial spins [2^ and twist s l29l . The formation of Q-lpalls after inflation have been extensively 
investigated in both GRV-M models [30] and GAU-M models [31, 32], in which SUSY is broken by either 
gravity or gauge interactions. As we will show, the Q-ball formation involves a nonequilibrium dynamics 
which is related to reheating problem in cosmology. 

The reheating process after the inflation period involves nonlinear, out-of-equilibrium, and nonpertur- 
bative phenomena so that it is extremely hard to construct a theory for the whole process, see the 2 
particle irreducible effective action as a remarkable approach [33] . In the first stage of reheating [pre- 
heating) ^ it is currently well known that the fluctuations at low momenta are amplified, which leads to 
explosive particle production. After preheating, the subsequent stages towards equilibrium are described 
by the wave kinetic theory of turbulence; Micha et. al. |34| recently estimated the reheating time and 
temperaute. These turbulent regimes appear in a large variety of nonequilibrium process, and indeed, 
the evolution of Q-ball formation experiences the active turbulence at which stage, many bubbles collide 
as observed in the next stage of tachyonic preheating [35] . During this bubble-collision stage within the 
reheating scenariOjJt is believed that gravitational waves may be emitted from the stochastic motion of 
heavy objects [H, [13] ■ The problem of gravitational wave emissions has been discussed only in the frag- 
mentation stage of Q-ball formation so far [H] , but not in the collision stage as opposed to the preheating 
cases. 

In this paper, we show analytically and numerically that in GRV-M and GAU-M models the approxi- 
mate trajectory of the AD fields is, respectively, precessing spiral or shrinking trefoil due to the quantum, 
nonrenormalisable, and the Hubble expansion effects. Moreover, we explicitly present the exponential 
growth of the linear spatial perturbations in both models. By introducing 3 -f 1 (and 2 -f l)-dimensional 
lattice simulations, we identify that the evolution in Q-ball formation involves nonequilibrium dynamics, 
including turbulent stages. Following the pioneering work on the turbulent thcrmalisation by Micha 
et. al. |34l |. we obtain scaling laws for the evolution of variances during the Q-ball formation. 

The paper is divided as follows. We explore both analytically and numerically the dynamics of the 
AD field in Sec. [Ill In Sec. IIIIl we study the late evolution of the AD fields and the process of Q-ball 
formation, introducing detailed numerical lattice results. Finally, we conclude and discuss our results in 
Sec. IIVI Two appendices are included. We obtain the equations of motion and their perturbed equations 
for multiple scalar fields in Appendix 1X1 In Appendix [Bl we find elliptic forms for the orbits of AD fields. 

II. THE AFFLECK-DINE DYNAMICS 

In this section, we investigate the equation for the orbit of an AD condensate. This orbit coincides 
with the well-known orbit equation in the centre force problem in classical dynamics, i.e. planetary 
motions so that we call the AD condensate "AD planet" sometimes. For the bound orbits, the effective 
potential should satisfy the condition where the curvature at the minimum of the effective potentials 
should be positive. In the presence of Hubble expansion, the effective potential depends on time; thus, 
the full solution of the orbit equations can be obtained numerically except the case that the AD field 
is trapped by a quadratic potential. In Appendix [Bl we obtain the exact orbit in this exceptional case 
when Hubble expansion is assumed to be small and adiabatic. The orbit of the AD planet, or more 
precisely an eccentricity of the elliptic motion in the complex field-space, is determined by the initial 
charge and energy density. In order to obtain analytic expressions of the orbit in more general potential 
cases in which we are more interested, we restrict ourself to work in Minkowski spacetime and on the 
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orbit which should be nearly circular. In Appendix [B] we also obtain the perturbed orbit equation and 
necessary conditions for closed orbits where the orbits come back to their original positions after some 
rotations around the minimum of the effective potential. In Bertrand's theorem [39j |. there are only two 
potential forms allowed to be closed orbits: isotropic harmonic force and the inverse-squared force. Each 
of the central forces gives dynamical symmetries, namely Fradkin tensor [ioj and Runge-Lenz vector l4ll . 
respectively. These dynamical charges are obtained both classically by the algebra of Poisson bracket [42| 
and quantum- mechanically by the corresponding Lie algebra in the abelian case [43i | as well as nonabelian 
case [4J]. By approximating phenomenologically motivated models that appear in the MSSM and using 
the results in Appendix [Bl we present, in this section, analytic motions of the nearly circular orbits and 
the pressure of the AD planets. Further, we check these analytic results with full numerical solutions. 

Let us consider a motion of AD fields in an expanding universe with scale factor a{t) and Hubble 
parameter H = d/a, where a over-dot denotes the time derivative. We investigate the AD field after they 
start to rotate around the origin of the effective potentials and the value of the U(l) charge pg is fixed due 
to negligible contributions from A-terms. By decomposing the complex (AD) field (j) as 4>{t) = cr(t)e**'^*', 
where a and 9 are real scalar fields, the equations of motion for a{t) and 9{t) (see Eqs. (|A8[ IA9p in 
Appendix [X| are 



3H& 



da 



0, 

3H9+-&9 = ^ 



dt 



= 0, 



(1) 
(2) 



where the conserved comoving charge density is defined by pq ~ a^a^O, and the effective potentials are 

2 

V± = V{a) ± Note that we will use shortly. From Eq. (jA10|) . the energy density pE and 

pressure p are given by 



1 .2 



1 .2 

p — —a 
^ 2 



V- 



(3) 



With various values of the charge density pq , Fig. [T] shows typical effective potentials V+ in Minkowski 
spacetime where we set a = H = 1. The models shown in Fig. [1] will be used later. 
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FIG. 1: (Color online) We show the effective potentials, V+ = V{a) + against a in two types of potentials 
which we call gravity mediated model (GRV-M Model) on left and gauge mediated model (GAU-M Model) on 
right. The potential in GRV-M Model has the following form, V{a) = icr^ (l - \K\ \na^) + bla'^, where, we set 
\K\ = 0.1 and 6^ ^ M ^ 9.20 x 10"^ The potential in GAU-M Model is V{a) = In (l + ct^) + h^a^, where we 
set 6^ ~ lO"'^". We choose the following values of pq: red-solid line for pq ~ 2.36 x 10~^ and green-dashed line 
for pQ — 1/e ~ 3.68 x 10~^ in GRV-M Model and red-solid line for pq ~ L40 x 10^ and green-dashed line for 
PQ ~ 1.41 X 10^ in GAU-M Model. 

Given an initial charge and energy density (or equivalently initial momenta and positions), the AD 
field oscillates around the value Ucr, which is defined by 
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where the orbit becomes circular when it starts from there, i.e. cr(0) = acr, o'(O) = 0. This orbit is 
bounded when the curvature is positive 



> 0. (5) 



For example, given a power-law potential such that V = Ai cr' where Ai is a coupling constant and I is 
real and a power of the homogeneous field a, the condition given by Eq. ^ implies that bound orbits 
exist for / < —2, < Hf Ai > and for —2 < Z < if Ai < 0, where we used Eq. (j4|). Another example is 
the case that a scalar potential is logarithmic, i.e. y = A2 In cr where the coupling constant A2 is positive. 
In this case, Eq. ([5]) is automatically satisfied. We investigate these two cases in more detail in appendix 

m 

Let us rescale the field a{t) as a{t) ~ ( ) ^i^) where oq is the value of a{t) at an initial time. It 
follows that the equations of motion in Eqs. (HI [21) are 





3a^ 




(ao) 











^ dV{a) 



0, (6) 



da dt 

where we defined pq = a^B = a^^ pq, and the terms involving and a/a are negligible under the 
assumption of an adiabatic Hubble expansion, i.e. 1, a <^ a. 

By introducing a new variable, u{t) = l/a{t), and using the second expression in Eq. ([6]), the first 
expression in Eq. ([5]) becomes the well-known orbit equation in the centre force problem such that 

(Pu , 1 f a\^dV , 



Notice that J{u, t) depends on time caused by the Hubble expansion, whereas the time-dependence in J 



vanishes when the potential V is given by a quadratic mass term, iM^cr^, where M is a mass of the AD 



field. 



A. Model A and Model B for MSSM flat potentials 



Let us introduce two models which appear in the MSSM in which SUSY is broken due to either gravity 
or gauge interactions. The former case, so-called, the gravity-mediated (GRV-M) model, has a scalar 
potential 

V- = i„.v(l + ..-,.^)+-|,,", (8) 

where m is of order of the SUSY breaking scale, which could be the gravitino mass scale TO3/2 evaluated 
at the renormalisation scale [13. Also, A is a coupling constant for the nonrenormalisable term, which 
is suppressed by a high energy scale, e.g. the Planck scale rUpi ~ 10^^ GeV. Here, if is a factor from 
the gaugino-loop correction, whose value is typically K ~ —[0.01 — 0.1] when the flat direction does not 
have a large top quark component [1, S^; thus, we concentrate on the case oi K < from now on. The 
power n of the nonrenormalisable term depends on flat directions. As examles of the directions involving 
squarks, the u'^d'^d'^ direction has n = 6, whilst the u'^u'^d'^e'^ direction is n = 10. For \K\ <^ 0{1), the 

first two terms in Eq. ([H]) can be approximated by - then find that 

IVP A^ 

Via) c —a' + ^^a- for n>l (9) 

which we call 'Model A', where we set Af^ = m^M*'^' and M has a mass-dimension, ~ 1, since 

1 = 2 — 2\K\ for \K\ <C 0{\). For small values of cr, we confirm that the power I is not approximately 

2 — 2\K\, so we will find a value of I numerically in that case later. 

In another scenario in which SUSY is broken by gauge integration, so-called, gauge-mediated (GAU-M) 
model, the scalar potential has the curvature with the electroweak mass at a low energy scale, whilst it 
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grows logarithmically at the high energy scale (which means that the potential is nearly flat as similar 
as the case of / = in Eq. ([9])). Therefore, 

(10) 

where Ms is the messenger scale (~ 10"* GeV) above which the potential grows logarithmically and m0 
is the same scale as Mg- We, thus, set Mg = m,/, for later convenience. Then, the scalar potential at high 
energy scale is approximately given by psj 

V^ml\n(—)\^^a\ (11) 

What follows is that we assume the orbit of the AD condensate is determined by the high energy scale 
where o'er ^ calhng this case, Eq. (fTTjl . 'Model B'. 

Using the results in appendix iBl we obtain the following quantities, W, $ and (w) by assuming that 
the dominant contribution in Model A and B is, respectively, either power-law term or logarithmic term, 
each of which corresponds to the first term in Eqs. ((9| fTTj). respectively. Here, we defined $ as a phase 
difference when the radial field a goes from the minimum value through the maximum one and back to 
the same minimum point, see Eq. (|B30P : in addition, (w) is given by a value of the equation of state 
averaged over a rotation of the orbit, see Eq. (jB2ip . The sub-dominant terms (nonrenormalisable terms) 
perturb the orbits by introducing infinitesimally small quantities sa, b where the subscripts correspond 
to the names of models introduced above. Thus, the main contributions are either Eqs. (|B31[ [B32)) or 
Eqs. (lB34l[B35l) . 



1. Model A - V{a) = + -4^a" 

By recalling Eq. ([S]), we obtain the following relations for n > /: 

>.^' = '" + ^'f + (12) 

where we defined a positive parameter, = "[^!|^2^)'' jvf^n"'^ ^ which is assumed to be infinitesi- 
mally small. We also obtain /J^ ~ (? + 2) (l + ^^e^) > 0, where /3 is defined in Eq. (|B24p . Substituting 
(3 into Eqs. (|B30l |B2T|) . we obtain $ and (w): 

^ f I — n \ , ^ 

* - ^r=^ 1 + , (13) 



Vr+2 V 2(n + 2) 

il - 2) (l ..^5^) ^1^2 f^^ Aljn^l) ^ ^^^^ 



(/ + 2)(l + e^l) l + 2y n{n + 2)il 

From Eq. the orbits for / = 2 — 2\K\ ~ 2 are nearly closed, but it is perturbed by the nonrenormal- 
isable term involved with sa- It results in that the periapsis appears to precess where the precession rate 
can be obtained from Eq. (|12p . The reader should notice that the denominator of the term involving eA 
in the second expression of Eq. p4p has I — 2 ~ ~2|iir| ^ C(l), which implies that it would be possible 
to have the non- negligible contribution from the term, even though ^ C'(l)- From now on, we restrict 
ourself that this is not the case; therefore, the dominant contributions appear as the leading orders in 
Eqs. (iniini) and Eq. which correspond to Eqs. ((63111532)) and Eq. (|B33|l . From Eq. ^ with 

eA — 0, our results recover the result published in ^9j, (w) ~ 
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2. Model B - V{a) = ml In {cj/m^f + ^^a" 
By introducing another infinitesimally small positive parameter, e_B = ^^^^P^^a"^ <C 1, we find 

4> pi 

W^^^il + es), (15) 



1 - 2 In 1 




1 1 2(ri-2) 


1 + 2 In 1 







W) = , — T ;C —1- (16) 



Since we are working in the high-energy regime, Ucr 3> tti^, the pressure of the AD condensate is hkely 
to be negative, see Eq. p6)) . From the second expression for $ of Eq. (fTS]) . the orbits are not closed and 
it should look like trefoil, see Eq. (|B34p . 

In an expanding universe, the above orbits for Model A and B suffer from the Hubble damping so that 
the orbits are naively expected to be precessing spiral or shrinking trefoil in the field space, respectively. 

B. Numerical results 

In this subsection, we present numerical results to check the analytic results which we found in the 
previous subsection. To do so, we use the full potentials, Eqs. ([51 [TU|) . instead of Eqs. (^1 [TT|) . and then 
solve Eq. ID) numerically in Minkowski spacetime as well as in an expanding universe. We adopt the 
4th order Runge-Kutta method with various sets of initial conditions, such as pq and . Since our 
analytical work holds as long as ^ we are concern with the two cases: a nearly circular orbit 

with = 0.1 and a more elliptic orbit with = 0.3. First of all, we parametrise Eqs. ([51 [TU)) by 

introducing dimensionless variables: a = cr/Af*, = ^ — ^ = l^|e~^/4, i — mt, x = mx in GRV-M 

pi 

Model and a = a /Ms, 6^ ^^l^ , i = M^t, x = Af^x in GAU-M Model. Since we know that m - 10^ 

GeV, M* - 101° GeV, rupi -"lO^^ GeV; hence, we can set hi - 9.20 x lO'^ ~ 0(10-^) in GRV-M 
Model, where we choose \K\ = 0.1. Notice that these choices are same as the ones used in pjj. On 
the other hand, we know that ^ Ms ^ 10* GeV; hence, we can set ^ lO^'^" in GAU-M Model, 
where we choose A ^ 10~^ as used in the GRV-M case. Notice that we can obtain the rescaled charge 
density pq and energy density pE, such that pq — mM^pq, pE = rn^MlpE in GRV-M Model and 
Pq = M^pq, pe = M^pE in GAU-M Model. 

Therefore, our rescaled potentials in GRV-M and GAU-M models for a n = 6 flat-direction are, respec- 
tively, 

V = \a^{\~2\K\\na)^hlo^, (17) 

V = In (1 + 0-2) +52^6^ (^8) 

where we omit over-rings for simplicity. The reader should notice that these variables that appear within 
the rest of this sub-section are dimensionless. We can also obtain the ratio defined by an energy density 
divided by a mass multiplied by a charge density, where the mass corresponds to m or Mg in either 
GRV-M or GAU-M Model, respectively 

In order to obtain appropriate initial values of cr(0), (t(0) and 9(0) satisfying the conditions tA, cs <C 
0(1) and giving not too small charge densities, we shall show that we need to choose only the initial 
values of 9{{)) in both GRV-M and GAU-M models. First, by ignoring the nonrenormalisable term in 
Eq. ((ni) for GRV-M Model, we can obtain Ocr = exp {^-2[k\ (^^(0) + 1^1 ~ l)) := ^^(0) from Eq. (gl), 
where we set Ucr '■— "'(0), which implies that we are setting that the initial phase is 37r/2. Since a has 
the maximum value at CT = a^r, we can set cr(0) := £'^u(fi)\J e'^{0) - \K\/2 from Eq. ([Bl4)) . which imply 
that EA ~ l2bl(j^{G) from the definition. We notice that cr(0) > 0{l) for ^(0) < 0(1); hence, it breaks 
the condition, <C C>{1). We can also see that (t(0) ^ 0(1) for 0(0) 3> 0(1), so the charge density is 
suppressed exponentially. Therefore, we are concern with the following two cases: ^(0) = \/2 and 1.0, 
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which give, respectively, ca 1.20 x 10"", pq - 2.36 x lO^^ and ea ~ 1.58 x 10-^ pq - 3.68 x lO"!. 

Similarly, in GAU-M Model, we choose that cr„. = \[^^~^ ^(0)' ^(0) e^^^l - |^^(0) and es = 

126^(7^(0) from the definition of es. Here, we also set the initial phase is 37r/2 due to (Jcr ■= o'(O). With 
this fact and the approximation, cTcr ^ C>{1), we need to have ^(0) ^ 0{1). In addition, we should have 
cr(0) < O(IO^) due to the condition, es < C(l). Therefore, we choose ^(0) = \/2 x 10"^ and V2 x lO'^ 
which gives, respectively, es ~ 1-16 x 10~^^, pq ~ 1.40 x 10^ and es ^ 1-20 x 10"^'', pq ^ 1.41 x 10^. 

Upon the above initial conditions, we initiate the numerical simulations with 8 different sets of the 
initial values in GRV-M Model and GAU-M Model summarised in TABLE |T] where we call each of the 
parameter-sets 'SET-1, SET-2,..., and SET-8'. In Fig. [Tl we also show, with the various charges which we 
introduced above, the effective potentials V+ for the GRV-M potential given by Eq. ([TT]) in the left panel 
and for the GAU-M potential given by Eq. in the right panel. Our time-step, dt, in the numerical 
simulations is dt = 1.0 x lO""* in the GRV-M case and dt = 1.0 x 10"^ in the GAU-M case. 



SET 


Model 


6(0) 


a(0) 


PQ 


eA or es 




Pe/pq 


1 














0.1 


1.46 


2 




V2 


~ 4.09 X 10"^ 


~ 2.36 X IQ-^ 


~ 1.20 X 10" 


11 


0.3 


1.51 


3 


GRV-M 












0.1 


1.06 


4 




1.0 


~ 6.07 X 10"^ 


~ 3.68 X 10^^ 


~ 1.58 X 10 


-2 


0.3 


1.09 


5 














0.1 


4.00 X 10"' 


6 




V2 X 10-^ 


~ 9.95 


~ 1.40 X 10^ 


~ 1.16 X 10" 


23 


0.3 


4.03 X 10~^ 
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GAU-M 












0.1 


7.22 X 10-2 


8 




V2 X 10"^ 


~ 1.00 X 10^ 


~ 1.41 X 10^ 


~ 1.20 X 10" 


17 


0.3 


7.25 X 10"^ 



TABLE I: We show 8 different parameter sets in both the GRV-M and GAU-M cases, where we call each of the 
parameter-sets 'SET-1, SET-2,..., and SET-8'. The initial parameters of cr(0) and (t(0) can be obtained by the 
values of ^(O). We also set 0(0) = ^ in all cases, and show the values of for GRV-M Model and the values 
of es for GAU-M Model. By substituting these values and choosing the values of the third eccentricity — O.l 
and 0.3, we obtain the dimensionless energy-to-(mass multiplied by charge) ratios, Pe/pq- 



1. The orbit of an Affleck-Dine "planet" in Minkowski spacetime 

First, we present numerical results in Minkowski spacetime in order to check our analytical results. 
We then give ansatze which are motivated by our analytic solutions in an expanding universe in the next 
sub-subsection. 

The motion of ct^ (i) In Fig. [21 we show the numerical solutions using the GRV-M potential with 
Eq. (fT7|) (left) and using the GAU-M potential with Eq. (fTS]) (right), and compare them with the corre- 
sponding analytic solutions which are given by Eq. (|B16[) . Using the initial values whose parameter sets 
can be seen in TABLE U we plot the numeric and analytic solutions in Fig. [2l In the top-left panel, the 
numerical plots (red-plus dots for SET-1 and blue-cross dots for SET-2) have the same amplitudes as 
the analytical ones (gree-dashed hne for SET-1 and purple-dotted-dashed fine for SET-2), we, however, 
can see the significant differences for the frequencies of each oscillation. We notice that these discrepan- 
cies come from the artifact of our choice as Z = 2 — 21^1 in Eq. ([9]), since the choice is not appropriate 
for a < 0(1), recalling cr(0) - 4.09 x lO"^ SET-1 and SET-2. Shortly, we wiU obtain numerically 
this power Z, and the obtained semi-analytic solutions match with the numerical ones. With SET-3 and 
SET-4, we can see that a{0) is not so small as opposed to the previous cases, i.e. a{0) ~ 6.07 x 10"^; 
thus, in the left-bottom panel of Fig. [2] we can see a nice agreement between numerical plots (red-plus 
dots for SET-3 and blue-cross dots for SET-4) and analytic plots (skyblue-dotted-dashed line for SET-3 
and black-dotted hne for SET-4). 

Similarly, we show the numerical and analytic plots for the GAU-M potential in the right-panels of 
Fig. [2] using the parameter-sets: for SET-5 and SET-6 in the right-top panel and for SET-7 and SET-8 in 



8 



the right-bottom panel. By changing the values of the third eccentricity (see TABLE the numerical 
plots deviate slightly from our analytic lines in the right-top and right-bottom panels of Fig. [2] as we can 
expect; in particular, we can see that our analytic values of both the frequencies and amplitudes of a^ify 
are larger than the numerical ones, and this difference can be significantly reduced when the orbits of the 
AD planets is nearly circular with = 0.1. 

As we have seen in the left-top panel of Fig. [2l our analytic value, I = 1.8, in Eq. ([9]) are not good 
enough to reproduce the numerical solutions since a{t) <C 0{1). Therefore, we set a trial function, 
/(ct) = iff" -I- 6^CT®, where a numerical value a is found by using the 'fit' command in the numerical 
software called 'gnuplot'. We find that a — 1.86002 := I is the best value of a, where we fitted this trial 
function f{a) onto the numerical full potential in Eq. (fT7|) for the range of u G [1.0 x 10~^ — 1.0 x 10~^], 
recalling a{Q) ^ 4.09 x lO"'^ in SET-1 and SET-2. Using this value of a as the value of / instead of 
I = 1.8, we plot the semi-analytic evolution for cr^(t) in Fig. [3] (green-dashed line for SET-1 and purple- 
dotted-dashed line for SET-2) against the corresponding numerical plots (red-plus dots for SET-1 and 
blue-cross dots for SET-2). Now, our semi-analytic solutions match with the numerical solutions. 
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FIG. 2: (Color online) Using the parameter sets summarised in TABLE HI we plot the numerical evolution for 
(T^(t) in both GRV-M Model (left) and GAU-M Model(right). In all of the panels except the case for the left-top 
panel, the numerical plots (red-plus dots and blue-cross dots) agree well with the corresponding analytic lines, 
which are obtained from Sec. Ill Al The disagreements between the numerical and analytic plots in the left-top 
panel come from the artifact that the analytical estimated value, I — 1.8, in Eq. (|9]). 



The average values of w{t) Using Eqs. p4l [T6| . we show both numerical values {wnum) and (semi- 
)analytical values (wana) of the averaged equation of state in TABLE [Til For all cases, the AD condensate 
has a negative pressure and one can say that the numerical values are of the same order as analytic values. 



The values of $ In TABLE IIIIl we show the numerical and (semi-)analytic values of <J> in both 
GRV-M Model and GAU-M Model, which are analytically obtained in Sec. Ill Al Our analytical values 
agree well with the numerical values. These values suggest that the orbits in GRV-M Model and GAU-m 
model are nearly either elliptical or trefoil, respectively. 
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FIG. 3: (Color online) Substituting tlie numerical value, I = 1.86002, into Eq. ((9)l, we plot the semi-analytic 
evolution for (J^{t). Our semi-analytic solutions agree with the numerical solutions. 





GRV-M Model v.s. Model A 


GAU-M Model v.s. Model B 




SET-l|SET-2 


SET-3 


SET-4 


SET-5 


SET-6 


SET-7|SET-8 




-2.42 X 10"^ 


-4.47 X IQ-^ 


-4.45 X lO""" 


-6.43 X 10"' 


-6.45 X 10"' 


-8.00 X 10"' 


{Wana) 


-3.63 X 10^^ 


-5.00 X 10"^ 


-6.43 X 10"' 


-8.04 X 10"' 



TABLE II: Using Eqs. p4l I16|l . we show the both numerical values {wnum) and analytical values (wana) for 
the averaged equations of state. The values of (wana) in SET-1 and SET-2 are semi-analytically obtained by 
substituting I — 1.86002 into Eq. Q. For all cases, the AD condensate has a negative pressure, and these 
analytical estimates are reasonable against the numerical values. 



$ 


GRV-M Model v.s. Model A 


GAU-M Model v.s. Model B 




SET-1 


SET-2 


SET-3 


SET-4 


SET-5 


SET-6 


SET-7 


SET-8 




1.591 


1.590 


1.605 


1.604 


2.210 


2.206 


2.221 


2.217 


^ana 


1.612 (analytic) or 1.599 (semi-analytic) 


1.605 


2.221 


2.221 



TABLE III: We show the numerical and (semi-)analytic values of $ in both GRV-M Model and GAU-M Model, 
which are analytically obtained in Section [II Al 



2. The orbit of an Affleck-Dine "planet" in an expanding universe 

We carry out our numerical simulation in an expanding universe when the inflaton field, which is 
trapped by a quadratic potential, starts to coherently oscillate around the vacuum during the reheating 
era. Then the evolution of Hubble expansion, H{t), and scale factor, a{t), follow as an ordinary nonrel- 

/ n2/3 

ativistic (zero-pressure) matter, see Eq. (IB33[) . For / = 2, we find H = ^^^^^^^ and a{t) — ( ^7^] , 

where ao is given by the value of a(t) at t — and we set uq = 0.1. We also set the initial time as 
to = 4x 10^ for GRV-M Model and io = 4 x 10^ for GAU-M Model. Notice that with this choice of to our 
simulation starts from the same physical time because we rescaled the time by either m ^ 10^ GeV or 
Ms ~ 10^ GeV, respectively. We again solve the equation of motion, Eq. ([T|), numerically using the 4th 
order Runge-Kutta method and compare them with the following ansatze. In order to see the significant 
effects from the Hubble expansion, we use SET-3 in GRV-M Model and SET-7 in GAU-M Model as the 
initial parameters. 

In an expanding spacetime, one can guess that our analytical results in Minkowski spacetime should 
be changed. In particular, the amplitude of a{t) may decrease due to the Hubble damping as we saw in 




numerical plots with SET-1 + 

semi-analytic line with SET-1 

numerical plots with SET-2 -"^ GRV-M Model 
semi-analytic line with SET-2 
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the quadratic case in Appendix IB 11 and similarly the frequency W in Eq. ([5]) should be changed. Hence, 
the orbit of the AD planet is processing spiral or shrinking trefoil in either GRV-M or GAU-M Model as 
one can see jlB]- Let us give an ansatz for <T^{t), 

Here, we use the Minkowskian values of a and W, and will obtain the possible values of ai, 2 in both 
models. From Eqs. (HUHl), by ignoring the nonrenormalisable term and recalling a{t) — uq i^-j^j , we 
can find the following proportion relations: (Tcr{t) oc (t + tp)^'*^^'^^'' ~ (t + to)^^^*^^"'^'-* and W{t) oc 
it+to)~^^^ ^ (t+tQ)^^ in Model A, where we used Z = 2-2\K\. In Model B, we obtain dcr oc (t+<o)"^ 
and W{t) cx + to)"^. Therefore, we set ai = 2^\k\ ' ^'^ ~ ~ '^2 \ k \ Model A, and a\ =4, a-i = —2 
in Model B. We believe that our ansatze are valid as long as the nonrenormalisable term does not play a 
role and the frequency of the coherent rotation, 0{W{t)), is rapid compared to Hubble expansion rate, 
0{H). The latter restriction implies that the rotation time scale is much shorter than the time scale of 
the Hubble expansion, i.e. W^'^lt) > [13]. 

The motion of tT^(t) In Fig.lH we plot the evolution of a'^{t) with the numerical data (red-plus 
dots) for GRV-M Model (left) and for GAU-M Model (right) and with the analytic data (green-dotted 
lines) using our ansatze Eq. The readers should compare the Minkowskian cases of SET-3 (left 

bottom panel) and SET-7 (right bottom panel) in Fig. [2] with the corresponding expanding background 
cases. For both cases, the amplitudes of a'^{t) decrease in time as we expected, and our analytic plots 
excellently agree with the corresponding numerical results. In the left panel of Fig. [H the difference 
between the analytic line and the numeric plots arises for the late time. We believe that this comes from 
the artifact of the approximation on / = 2 — 2|iir| in GRV-M Model, Eq. p7|) . since the values of (7^{t) 
decrease to the region where the above approximation does not hold, i.e. for a <C C(l) as we saw in the 
left-top panel of Fig. 
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FIG. 4: (Color online) We plot the evolution of (T'^(t) with the numerical data (red-plus dots) for GRV-M 
Model (left) and for GAU-M Model (right) and with the analytic data (green-dotted lines) by using our ansatze 
introduced in Eq. (|19p . 



The motion of the equation of state: w{t) = p{t)/pE In Fig. El we plot the numerical values of 
the equation of state, which is given by w(t) = p{t)/pE, where p{t) and pe in Eq. ^ are the pressure 
and energy density of the AD condensate. The averaged pressure over the rotations seems to be negative 
in GRV-M Model, see the left panel; whereas, the pressure in GAU-M Model is always negative, see 
the right panel. The frequencies of the rotation for w{t) in both cases are, respectively, similar as the 
corresponding frequencies of a'^{t), see Fig. |4l however, the phases are different from the phases of a'^{t) 
by TT. 

In summary, we have analytically obtained the nearly circular orbits for both GRV-M Model and GAU- 
M Model in Eqs. ^T7\ [T5)) approximated by Eqs. ([51 [TT|l. We then checked that the (semi-)analytic results 
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FIG. 5: (Color online) Using the initial conditions called SET-3 (right-panel) and SET-7 (left-panel) in TABLE 
U we plot the numerical values of the equation of state which are given by w{t) = p{t)/pE, where p{t) and ps are 
the pressure and energy density of the AD condensate. 



in Eqs. [0)1 and Eqs. ^T5[ [TB)) and our ansatze in Eq. agree well with the corresponding numerical 
results obtained by numerically solving Eqs. pi IB12p . In the rest of this paper, we investigate the late 
evolution for the AD condensates once the spatial perturbations generated by quantum fluctuations or 
thermal noise from the early oscillation [lll | become non-negligible due to the negative pressure presented 
in TABLE n] and Fig. [5] 



III. Q-BALL FORMATION AND THE THERMALISATION IN MINKOWSKI SPACETIME 

In this section, we analyse the late evolution of the AD condensates in both GRV-M and GAU-M 
models, in which we find that the spatial perturbations are amplified exponentially due to the presence 
of the negative pressure, and the presence of negative pressure supports the existence of nontopological 
solitons, i.e. Q-balls. As a process of reheating the Universe, the dynamics of the Q-ball formation 
is nonequilibrium, nonperturbative, and nonlinear process, and it includes three distinct stages: pre- 
thermalisation (linear perturbation), driven turbulence (bubble collisions), and thermalisation towards 
thermal equilibrium. As opposed to the reheating process, we report that the driven turbulence stage 
lasts longer and the subsequent thermalisation process is different, which is caused by the presence of 
nontopological soliton solutions. During the turbulent stages, we find scaling laws for the variances of 
fields and for the spectra of the charge density. In addition, we adopt numerical lattice simulations to 
solve classical equations of motion, where our numerical code is developed from LATfield j48|], and we 
present detailed nonlinear and nonequilibrium dynamics (some videos are available [igj). 



A. Linear evolution — Pre-thermalisation 



The late evolution, after the AD condensate forms, depends on the properties of models. In the standard 
AD baryogenesis scenario the condensate govern by the quadratic potential, Eq. (|B1|) . decays into 
thermal plasma which may give our present baryons/lcptons in the Universe. By including quantum 
and/or thermal corrections in the mass term as in Eqs. ([51 llOp . the subsequent evolution may be different 
from the standard AD scenario since the AD condensate has a negative pressure. The negative pressure, 
which causes the attraction force among particles in the condensate, amplifies exponentially the linear 
spatial fluctuations. We see this exponential growth for the linear perturbations in nearly circular orbit 
cases with the growth rate Sm, and obtain the most amplified wave- number fc^, which give a rough 
estimate on the nonlinear time Iml and the radii of bubbles created just after the system enters into a 
nonlinear regime. As long as the perturbations are much smaller than the background field values, we 
call this initial linear perturbation stage ^ pre-thermalisation\ 
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1. Arbitrary and circular orbits 



Let us consider the linear spatial instability for an AD condensate in Minkowski spacetime. First, we 
perturb the AD field (j) with the linear fluctuations, Sa and SO. Equations of motion for da and 69 are 
given by Eqs. (|XTT1 ET2|, 

Sa - (v^ +9^ - V"^ Sa - 2a9S9 = 0, (20) 
SO + —69 - V^60 + % (aSa - aSa) = 0. (21) 



a 

Let us rescale Sa and SO in the following form 

6a - <5f7oe^'*^+^'''^, 60 - (56loe'5(*)+»k.x^ ^22) 

Notice that both of the exponents S{t) should be the same in terms of a function of the wave number k, 
because we are concerning with linear perturbations. Substituting Eq. ([22|) into Eqs. ((20l [2T|) . we obtain 



5*2 + k2 - + V" -20S 

2d(s-^] 52 + ^ + k2 




(23) 

where V" = and we ignore the terms S, assuming that the linear evolution is adiabatic, i.e. ^ S 
(WKB approximation). Notice that this assumption is violated only at the beginning of this linear 
evolution as we will see in the numerical subsection. Sec. IIII CI The nontrivial solution for S can be 
obtained by taking the determinant of the matrix in Eq. (|23p . namely 



F{S{k),k^) EE 54 + ^^3^(2k2 + 302 + y")^2 

fk2 - 302 + V") 5 + k2 (k2 -O^ + V") = 0. (24) 



a V 



Notice that the terms involving a vanish if the orbit of the AD field is exactly circular. By looking for 
the most amplified mode fc^, which is defined by ^p-|j.2 ~ from Eq. ([M]) . it implies that 

kfn = - 5 (^^ + ^) > 0, (25) 

where the inequality comes from the reality condition for fc,,„. By concerning with this mode in Eq. (|25p 
and by solving F{S{k),k^) in Eq. the solution for S„i = S{k = fc„i) is 



^ (5^2 - V") ± 29 J (92 - V"Y + 2 {^f (3^2 - V" 

Sm = — (26) 

2 4.9^-^) 



in which we are interested in the growing mode, i.e. Re{Sm) > 0. Substituting Eq. (p6| into Eq. (l25|l . we 
may obtain the most amplified mode. Although it is rather hard to analytically solve Eq. we know 
that the only one instability band exists for exactly circular orbits where a = 0; 

<-k^ <d^ -V"{a), (27) 

where and a = a^r are time-independent due to the circular orbits. 

In addition, we can estimate a possible nonlinear time ijvi when the spatial averaged variance, Yar{a), 
becomes comparably large to the corresponding homogeneous mode a. Here, we defined Var (cr) = 

((7(x, t) — ct)^, and a hat and a bar denote an original field and a spatial average of the field, respectively. 
Notice that the nonlinear time in [l^, is defined by the time when the linear fiuctuation 6a for the 
most amplified mode becomes comparable large to the homogeneous-mode; however, our definition is of 
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more advantages as we see in the numerical subsection, Sec. IIII CI The nonlinear time with our definition 
can be given by 



Var{a) ^ Sa^ exp (^2iv(^)T^ ^ ^^o exp (^^ 2 ^5„.)^ - al (28) 



O tNL ~ ^, + In ( ^ ) . (29) 



Here, we approximated that ^ i^Smj and that the orbits over rotations with the period r, 

Eq. (jB17[) . can be expressed by the integral form as shown in Eq. (pS)) . As we assumed, the spatial 
averaged variance of this field is not fully developed over all modes except k = km until t ^ t^,, where 

is a typical time scale when the variance starts to grow with the growth rate ^S*. 

Our main interest in this pre-thermalisation stage is the evolution of the number of particles in terms 
of modes, so that we consider pg as the particle number here. For a free field theory, both of the positive 
and charged particle occupation number develop equally. The present case, however, gives different 
consequences due to the presences of nonlinear interactions and the initial inequality of a charge density 
(baryon asymmetry). Without loss of generality, we can focus on the case where the positive charge 

is initially present. Since the charge density is given by pg ~ a^O, we can approximately obtain the 
evolution in the linear regime using Eqs. HUll), 

PQ~a^{t)V^Se. (30) 

Hence, the charge density evolves due to the linear fluctuation of the phase field. Let nf{t) be the 
amplitude of Fourier-transformed positive and negative charge density, (x, t) , which are defined through 



k ' 



the following decomposition, pg ~ ?^ (x, t) — n (x, t). Notice that the Fourier transformed functions, 
are related to, but are potentially different from the corresponding quantum mechanical expressions, nj, = 
a^ak, = b^^bk and Q — J (Pxpg — J (J'^yl'^ i^k ~ "fe )■ Here, are occupation numbers for positive 

and negative charged particles in a free field theory, and a^, a|,, bk and b^. are the annihilation/creation 
operators for both of the particles, respectively. Since we are interested in the growing mode for the 
positive charge density n'^{t) in Eq. (|30p which is initially zero except the zero- momentum mode, it 
implies that using Eq. ((22)) 

n+{t) ~ k^\6e„\ /'dtV(f)e(^W>*, 



Jtn 

e\Seo\al / , (xe(^) (31) 

where to is a numerical value and we assumed cr^(i) ~ cr'^j., going from the first line to the second 
one. Therefore, the evolution of the positive charged particle number for a mode k is proportional to 

g(S(/c))(t-to)_ 

Our results, Eqs. ([25l [26)) . are generalisations of the known results [s^, 'HT'l, in which the orbit of the 
Affleck-Dine condensate is exactly circular. We also obtained the nonlinear time Inl in Eq. (|^5)) and the 
exponential growth of the particle number in Eq. pip . 



2. Nearly circular orbits in Model A and B 

Using the results obtained in the previous subsection, we can compute the most amplified mode (k^j) 
and the growing mode/s'm^ averaged over one rotation of the nearly circular orbits for the models 

introduced in Section fll Al i.e. Model A and Model B. We shall confirm that these values are the same 
as the cases when the orbits are exactly circular, which implies that the instability band, Eq. (|27)) . could 
exist even for the present nearly circular orbit cases. 



14 



Model A: Substituting the expressions, ct/ct, 9"^ and V" [c.f. Eqs. (|BT61[BT8)) and Eq. ((HI)), into 
Eq. (|26l). we obtain the averaged growing factor for Model A where IvP > 0: 



^"Z " ^^^V ^ + 2(n + 2)a-2) ^-^j ' (^2) 
(k„) ^ (^1 + (^^2)(/-2)(Z + 6) '^j ' 

where we substituted Eq. ([32]) into Eq. ([25| to obtain (k^) and these results are consistent with the case 
for the exactly circular orbit. In order to satisfy (k^„) > 0, we should have / < —6, < ^ < 2, and 
Eq. (j32p implies that the condensate is unstable against spatial fluctuations when the pressure is negative 
with < ; < 2, see Eq. (|B32)) . 

We can check the results ^\ that (s.^^^ ~ ^ [l + and (k2„) ~ m'^\K\ {l ~ by setting 



I — 2 — 2\K\ in Eqs. ((32l [33)1 and by ignoring the nonrenormalisable term as did in [5l|, i.e. {^Sr, 
■^~2~~ ^) ^c^'^' ^i^d (k^j) ~ |iir|A/^ ^1 — -^^^ CTcr^'^'. These are of the same order as their results, 
recalling that aZr^^^ - 0(1) due to \K\ < 0(1). 

Model B: Similarly, we can also obtain the averaged growing factor for Model B from Eq. (fT5|) 



" 71^ I' " j ' ^^"^ 2^ I' " 3(^^^ j 

whose leading orders reproduce the results [s^ in which case that the AD orbit is assumed to be exactly 
circular and ignored the nonrenormalisable term. 

Before we finish this subsection, let us remark the classical and absolutely stability of AD condensates. 
Lee found [l3l that the dispersion relations for the waves of linear fluctuations from Eq. ((24|) when the 
orbits of the AD fleld are bounded. For the longwave-length limits, there exist one massive and one 
massless modes. The massless mode can be interpreted as the sound wave whose sound speed should be 
real for the classical stability reason, and the squared value of the sound speed is related to the value 
of {w) in Eq. (|B21|) . Therefore, this stability condition for the sound waves corresponds to the sign of 
the pressure in the AD condensates. In other words, the AD condensate has a negative pressure when 
the sound speed is imaginary; equivalently, it is classically unstable against spatial fluctuations. The 
zero-pressure AD condensate whose energy density is minimised with respect to any degrees of freedom 
is equivalent to the Q-matter phase as Coleman discussed in [l^, where the absolutely stable Q- matter 
can be excited by classically stable sound waves. 



B. Non-linear evolution and nonequilibrium dynamics 



1. Driven (Stationary) and free turbulence 



Even when the perturbations are fully developed to support the nonlinear solutions, the system is 
still far from equilibrium. Indeed, the system enters into more stochastic stages, 'turbulence regimes', 
where the strength of the turbulent behaviour depends on the "Reynolds" number [s^l- As a theory of 
reheating Universe, a general nonequilibrium system goes through two different turbulence stages, going 
from driven turbulence to free turbulence stages. A major energy transfer from the zero mode takes 
place during driven turbulence. Garcia-Bellido et. al. [37| observed that bubbles form and collide during 
this stage in tachyonic preheating, and they suggested that the bubble collisions can be an active source 
of gravitational waves [53|. In the usual reheating scenarios, this stage terminates when the energy left 
out in the zero-mode becomes smaller than the energy stored in other modes (created particles). Since 
the energy exchange between zero-mode and other modes becomes negligible, the particle distribution 
is self-similar in time (free turbulence) and evolves towards thermal equilibrium. In the free turbulence 
stage, the quantum effects change the late evolution significantly, and the created particles are distributed 
followed by Bose-Einstein statistics rather than by a classical manner. As long as an active and stable 
energy source exists in momentum space, we expect that the driven turbulence stage lasts for a long 
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time. In the case of Q-ball formation, we expect that this active energy source corresponds to the excited 
states of Q-balls; hence, we can expect that the driven turbulence stage lasts longer compared to the 
linear perturbation regime as opposed to the usual reheating Universe scenarios. Note that during this 
thermalisation stage the transition from the classical to quantum regime becomes important [s^l ; in the 
rest of this paper we concentrate on the case where the system is govern by classical evolution all the 



time. In turbulent stages, the scaling law can be found 3J|: 



Var(cr) cx t^, (35) 

where the power p depends on the parameters of models, e.g. the relativistic values of p are p = 2m-\ 
in the driven turbulence regime and p — — 2m-\ ™ '^'^ ^'^^^ turbulence regime. Here, m is the number 
with which particles interact. For the free turbulence regime, the particle number distribution follows a 
scaling law from the time i/r-ee when the free turbulence turns on, namely 

4 

"-fc(i) = = t/ree), (36) 

where fc, = kt~ . 



Thermal equilibrium state in the presence of nontopological solitons 

In this sub-subsection, we show that the condition of the negative pressure is the same as the existence 
condition of Q-balls, known in p^ . This does not always mean that the spatially unstable condensate 
evolves towards to Q-balls; with given initial conditions, the condensate may evolve into other thermo- 
dynamically favour states in which the free energy is minimised. 

The ansatz of Q-balls claims that 9, which corresponds to the "chemical potential" cj, is constant, 
and that the radial field a should be time-independent and depend on the radius r of the Q-ball, i.e. 
(j) = (T(r)e"^*. Hence, the existence condition of Q-balls at zero-temperature is 



(37) 



This condition implies that the potential should grow less quickly than a quadratic term; thus, it is 
equivalent to the fact that the AD condensate has a negative pressure for Z < 2 in Eq. ^ , see Eq. (jB32|) . 
Notice that this condition only tells that Q-balls may appear after a unstable AD condensate fragments. 
The evolution to the thermal equilibrium state is rather hard to compute analytically, and it is related to 
stability problems of the Q-balls [l9|, jl^]. Therefore, we conduct numerical lattice simulations that give 
the entire processes of nonlinear as well as nonequilibrium evolution. 



C. Numerical results 



In this subsection, we present detailed numerical results involved with lattice simulations for both 
GRV-M and GAU-M models with the parameter sets, SET-3 and SET-7 shown in TABLE [H we then 
check our analytical results obtained in the previous sections. In order to solve the second-order partial 

differential equations, ^ — V^<^ -I- ^ = 0, with the potentials introduced in Eqs. (flTl fT8|) . we use the 
following appropriate parameters: dx = 0.2, dt = 0.02 in GRV-M Model and dx = 5.0, dt = 0.2 in 
GAU-M Model, which minimise the numerical errors. Here, dx is the fundamental lattice space and dt 
is the time step. Note that the variables in this subsection are normalized by appropriate energy scales 
as in Sec. IIIBI We then conduct 3-1-1 (and 2 -|- l)-dimensional lattice simulations with 512"^ (and 512^) 
lattice units, imposing a periodic boundary condition. Our initial conditions are, 4>o — (t^o + (500 and 

00 = 00 + (500 J where the initial fluctuations, (50o and 50o, are of a Gaussian noise, which are responsible 
for "quantum" fluctuations. Their fluctuations, S(po and (50o, are of order 10~^ in GRV-M case and of 
order 10^'^ in GAU-M case. In order to visuahse these detailed evolution, we use a 3D software, 'VAPOR' 
[s^ . and some videos of our numerical results are available in (49j . 
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1. Pre-thermalisation 

The initial evolution — Non-adiabaticity: In the top two panels of Fig. [6l we plot the amplitude 
of [t) , where we took the average of (t) over the axes of k. We show the amplitudes of (t) for 
GRV-M Model in the left panel and for GAU-M Model in the right panel with two different time steps. In 
the panels, we indicate the analytical values of the most amplified modes fc,„ obtained from Eqs. ([33l [34| 
with black-dashed vertical lines. In GRV-M Model, the amplitude with t — 30 (green-dashed line) is a 
little noisy to see the first peak ki in terms of k. Our analytical estimate, km ~ 2.88 x 10^^, is located 
at a more infrared region than the point k — ki ^ 0.34, and the periodic structure can be seen in the 
higher-momentum space. In GAU-M Model, on the other hand, we can confirm that our analytical value, 
k = km ^ 1-22 X 10~^, agrees with the numerical value, fci ~ 1.70 x 10~^, in the green-dashed line; 
however, the analytical value appears in a slightly more infrared region. We also observe the periodic 
structure in the higher- momentum modes as it was reported in [34I. In the middle panels (GRV-M 
Model on left and GAU-M Model on right), we compare both the zero- mode, cr^ (red-solid hues), and the 
homogeneous field, cr^ (green-plus dots), shown in the bottom panels of Fig. O The middle panels in both 
cases show that the zero-mode does not decay quickly, and it oscillates around cr^ — tr^^. We can also 
check that our numerical parameters are appropriate, minimising numerical errors. In the bottom panels 
of Fig. El we plot the evolution of nj, {t) for the modes both km (red-solid lines) and ki (green-dashed 
lines). In the left bottom panel, we can see the exponential growth of the amplitude in GRV-M Model for 
both modes, and step-like particle production exists at the beginning of the evolution as broad resonant 
preheating f3| (c.f. Eq. ([30]) ). and it begins to create the particles exponentially afterwards. The particles 
are produced quickly when the zero- mode cr^(i) increases in time at the beginning, see the middle panels. 
This is the different feature of the evolution compared to the case of resonant preheating, where particle 
production for the broad resonance occurs nonadiabatically when the zero mode (inflaton field) crosses 
the zero axis. In the right bottom panel, we can see more clearly the step-like particle creations for both 
modes, and then this step-like evolution smoothes out, which leads to the exponential particle production 
as in the GRV-M case. We believe that the adiabatic condition, S <^ S^, is "softly" violated only in this 
initial stage since we can not see the clear exponential growth at the beginning of this evolution. In the 
next paragraph, we discuss the late linear evolution when this nonadiabatic evolution ceases, and show 
that our analytical results agree with numerical ones more nicely. 

Up to the nonlinear time: In Fig. [71 we show the evolution of the various physical quantities in 

the late stage of linear perturbations: n^, tr^ and Var(cr). The top panels plot the amplitude of with 
various times in both GRV-M Model (left) and GAU-M Model (right). Notice that we plot it against the 
logarithmic scale of k as opposed to the linear scale seen in the top panels of Fig. [51 For all time steps 
shown there, our analytical values of km agree well with the first peak mode ki, at which the amplitudes 
are most amplified. Notice that the zero-momentum mode does not decay in both cases. After the first 
peak of the amplitude is well developed, the second peak appears in the spectra, and later the third 
peak can be barely observed. Roughly speaking, the n*'' peaks appear around the values which are km 
multiplied by n. These higher peaks are suppressed by rescattering processes in which a particle from the 
first peak transfers some of its momentum to a particle from the zero-momentum modes (AD condensates) 
[ssl ]. Later, all modes of the particle spectra, n^, develop quickly, but the first peak is still visible. The 
middle panels illustrate the evolution of a zero-mode field and the variance of the field Var{a) up to 
the nonlinear time t = ^atl- As we saw in the top panels, the zero mode does not decay even after the 
nonlinearity comes in, whilst the variance of the field develops exponentially from t ~ 140 in GRV-M 
Model (left) and from t ~ 600 in GAU-M Model (right). This delay of the exponential growth comes from 
the fact that the other modes do not evolve initially except the mode km', thus, we can set these times 
as defined in Eq. (j^ . We fit a function, cx exp (2Snum{t — t^)j , against the exponential evolution 

for the variations, where Snum is a numerical value, and we obtained Snum ~ 4.45 x 10^^ in GRV-M 
Model and Snum ^ 6.72 x 10~^ in GAU-M Model, which match satisfactorily with the analytical ones in 

Eqs. (121 [M]), where we computed as (Sm"^ ~ 4.20 x 10"^ in GRV-M Model and (Sm"^ ~ 7.07 x lO'^ 

in GAU-M Model. From the middle panels, the nonlinear time is approximately both t^L ~ 420 in 
GRV-M Model and t^L ^ 2200 in GAU-M Model, and these values agree well with the analytical 
estimates in Eq. (pS)) . where the analytical values are ^atl ~ 262 -I- 140 '-^ 422 in GRV-M Model and 
tNL ~ 1628 -I- 600 ~ 2228. In the bottom panels, we plot the evolution of the amplitude for the 
first peak mode (red-plus dots), second peak mode (green-cross dots) and the analytical most amplified 
modes (purple squared-cross dots). The numerical values of the exponents for the most amplified modes 
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FIG. 6: (Color online) In the top two panels, we plot the amplitude of ii^{t) with two different time steps for 
GRV-M Model in the left panel and GAU-M Model in the right panel, where we took the average of nj^ {t) over the 
axes of k. The black-dashed vertical lines indicate the analytical values of the most amplified modes km obtained 
from Eqs. (|33l|3U. In the middle panels (GRV-M Model on left and GAU-M Model on right), we compare the 
zero-mode (red-solid lines) and the homogeneous field (green-plus dots) obtained in the bottom panels of 
Fig. [2] In the bottom panels of Fig. (6] we plot the evolution of n^(i) for both analytic values km (red-solid lines) 
and numerical values fci (green-dashed lines) of shown in the top two panels. 



fc™ in blue long-dotted lines, {Snum ~ 4.55 x 10"^ in GRV-M Model and Snum ^ 7.11 x 10"^ in GAU-M 
Model) match with the analytical ones in Eqs. ([32l [34 | . The second peaks k2 in black short-dotted lines 
start to grow at t ~ 220 in GRV-M Model and at t ^ 1300 in GAU-M Model, and we can set these 
values as to defined in Eq. (|?1]) . The initial behaviour of the amplitude of second peak of seems to be 
quasi-periodic, which implies that (S) for the mode, A:2, is pure imaginary, see Eq. (|22p (c.f. the bottom 
panels of FIG. 5 in [sj). Surprisingly, the growth rates are about twice larger than the values of both 

Sm) and Snum for Var((T) and fci. Note that the initial evolution is not adiabatic, so that the growth 
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rates are not strictly exponential as we have seen in the bottom panels of Fig. [51 For example, the growth 
of the first peaks, km (or fci), in GAU-M Model is not exponential initially, but it becomes exponential 
as the growth of the second peak mode k2. 
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FIG. 7: (Color online) The top panels plot the amplitude of with various times in both GRV-M Model (left) 
and GAU-M Model (right). The analytical value of the most amplified mode km in black-dashed vertical lines 
agrees with the first peak, fei, of the spectra in both cases. The middle panels show the evolution of zero-mode 
field, (T^ (red-plus dots), and the variance of the field, Var((T) (green-cross dots), up to the nonlinear time t — Inl, 
where we can set Inl ~ 420 in GRV-M Model and Inl ^ 2200 in GAU-M Model. In the bottom panels, we plot 
the evolution of the amplitude for the first fci (red-plus dots), second peak k2 (green-cross dots) modes and 
the analytical most amplified modes km (purple squared-cross dots). 



Bubbles pinched out of filaments: In Fig. [51 we show the snapshots of the positive charge density 
n+(x) for GRV-M Model (left panels) and GAU-M Model (right panels) around t^tML, where 'Timestep' 
in the panels denotes the actual time divided by 10 in GRV-M Model and the actual time divided by 
10^ in GAU-M Model. The colour bars illustrate the values of the positive charge density. We can see 
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long-wavelength objects (sometimes called 'filaments') in both cases, and the charge in some regions is 
compactified into spheres, sec bottom panels. These filaments and bubbles correspond to the nonlinear 
solutions, which may be nontopological strings [56| and the excited states of Q-balls, respectively. The 
radii of these bubbles are of the same order as the wave-length which corresponds to the most amplified 
modes, km- As we will see in the next subsection, these bubbles grow by colliding and merging each other. 
Note that this bubble creation is nothing to do with bubble nucleation in first-order phase transition 
as opposed to the case in [iJl, in which case the AD condensate is classically stable against special 
perturbation, but not quantum mechanically. 




FIG. 8: (Color online) In the top and bottom panels, we show the snapshots of the positive charge density 
n^(x) for GRV-M Model (left panels) and GAU-M Model (right panels) around t ~ ijvL, where 'Timestep' in the 
panels denotes the actual time divided by 10 in GRV-M Model and the actual time divided by 10^ in GAU-M 
Model, and the colour bars illustrate the values of the positive charge density. After the nonlinearity is fully 
developed, many bubbles form, which are pinched out of "highly" charged filaments. 



2. Nonlinear evolution 

Bubble collisions and mergers: In Fig. [9l we show the snapshots of the positive charge density 
for GRV-M Model in different time steps up to i = 6000, where 'Timestep' in the figure denotes the 
actual simulation time divided by 10^ and the colour bars illustrate the values of the positive charge 
density. After the system goes into a nonlinear regime, we can see a few lumps in the first few panels of 
the snapshots, and those lumps merge into larger lumpy objects. Finally, we can see a large cluster which 
consists of a complicated inner structure, see the last snapshot. Recall that we are using the periodic 
boundary condition. 

Fig. [TUl shows the detailed evolution of the positive charge density for GAU-M Model in different time 
steps up to t = 60000, where 'Timestep' in the figure denotes the actual simulation time divided by 10^ 
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FIG. 9: (Color online) We show the snapshots of the positive charge density for GRV-M Model in different time 
steps {t = 1000, 2000, 3000, 4000, 5000 and 6000), where 'Timestep' in the figure denotes the actual simulation 
time divided by 10^ and the colour bars illustrate the values of the positive charge density. A few created lumps 
collide and merge into a large cluster in the end. 



and the colour bars illustrate the values of the positive charge density. A large number of small bubbles 
can be observed, and nearby bubbles collide and merge into larger bubbles. In the final panel, there are 
smaller number of bubbles left (compare to the first panel) . We believe that this time arrow is followed 
by the fact that the total energy of large bubbles is smaller than the total energy of smaller bubbles, c.f. 
fission stability of Q-balls in [l^ . These large bubbles are able to carry a large amount of charge inside 
of them as we saw in the left panel of FIG. 9 in [l^ in the "thin-wall" Q-ball limit. 

The differences of the evolution between GRV-M and GAU-M models come from a number of facts, 
e.g. different initial conditions, stability conditions and momentum fluxes due to asymptotic profiles at 
a large distance from the cores, see [l9| . 

Distributions of the negative charge density: We show the snapshots of the negative charge 
density for GRV-M Model (left panel) at t = 6000 and GAU-M Model (right panel) at t = 1.0 x 10^ in 
Fig. where the colour bars illustrate the values of the negative charge density. These times correspond 
to the same times as in the final snapshots of Figs. [51 UHl The values of charge density in both models are 
much smaller than the values of positive charge density in Figs. |9l 1101 This implies that we are observing 
the plots of thermal plasma rather than the charged (nonlinear) lumps. Their distributions are quite 
different each other. The negative charge density for GRV-M Model is surrounded by the large positive 
charged cluster seen in the last panel of Fig. [SI and it is distributed all over the lattice; whereas, for 
GAU-M Model the distributions of the negative charged plasma are highly concentrated only around the 
surface of the lumps (compare the last panel of Fig. [TO)) . 

Driven turbulence: The top panels of Fig. [T^show the evolution of zero-mode (red-solid lines) and 
the variations for a (dotted-dashed purple lines), whose latter evolution are fitted by a function, oc t^^, 
(black dashed lines) , where 71 is a numerical value as the power of Eq. ([55)1 . For both models (GRV-M 
Model on the left panel and GAU-M Model on the right panel), the asymptotic evolution after the linear 
perturbation regime is overlapped by the function, where 71 ~ 0.121 for GRV-M Model and 71 ^ 0.235 
for GAU-M Model. Our analytic values can be matched by setting p ~ 0.111 with m = 5 in GRV-M 
Model and p ^ 0.250 with to = 3 in GAU-M Model, see Eq. ([55]) . Hence, we could identify this regime as 
driven (stationary) turbulence, and the main dynamics in each model is caused by either a five-particle 
interaction or three particle interaction, respectively. Note that our nonrenormalisation term has a 0^ 
term in both models. In the middle and bottom panels of Fig. [121 we plot, respectively, the amplitudes of 
n'^ and n^T in different times for GRV-M Model (left panels) and GAU-M Model (right panels) . For of 
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FIG. 10: (Color online) We illustrate the detailed evolution of the positive charge density for GAU-M Model in 
different time steps [t = 10000, 20000, 30000, 40000, 50000 and 60000), where 'Timestep' in the figure denotes 
the actual simulation time divided by 10'^ and the colour bars illustrate the values of the positive charge density. 
There are smaller number of bubbles left in the end. 




FIG. 11; (Color online) We present the snapshots of the negative charge density for GRV-M Model (left 
panel) at t = 6.0 x 10'' and GAU-M Model (right panel) at t = 6.0 x lO'*, where the colour bars illustrate the 
values of the negative charge density. The negative charge density for GRV-M Model is surrounded by the large 
positive charged cluster; however, the distribution spreads out over the lattice space, whereas the negative charge 
density for GAU-M Model is concentrated around the positive charged lumps (compare them to the last panels 
of Figs.[9l[l0j. 

GRV-M Model, the amplitudes of the high momentum modes grow in time, whilst the lower momentum 
modes do not decay completely and stay for a long time. We fit a function, cx k^^^, (yellow dotted lines) 
where 72 is a numerical value onto the spectra at i = 6700 for the region where the function is fitted as 
shown in black dashed lines. We find that 72 ~ 1.62 for the n'^ case and 72 ~ 0.37 for the case. In 
the right middle and bottom panels, we plot the amplitudes of for GAU-M Model in various times. 
The amplitudes of the high momentum modes decrease as opposed to the GRV-M case, and the slopes 
of the spectra for nfa.tt = 63000 in yellow-dotted lines are steeper than the GRV-M case, where we fit 
the numerical spectra by the following values shown in black dashed lines: 72 ^ 3.95 for the case and 
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72 ~ 1-74 for the n^. case. 
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FIG. 12: (Color online) Left panels (GRV-M Model) and right panels (GAU-M Model): the top panels show 
the evolution of zero- mode (red-solid lines) and the variations for a (dotted-dashed purple lines), whose latter 
evolution are fitted by a function, oc t'''^ , (black dashed lines) where 71 is a numerical value as the power of 
Eq. (I35p . In the middle and bottom panels, we plot, respectively, the amplitudes of and in difi'erent times 
for both models, and we fit them by a function of oc k where 72 is a numerical value. 



3. From driven turbulence to near equilibrium - Thermalisation: 



In order to significantly reduce the simulation time, we carry out 2 + 1-dimensional lattice simulations 
with the same initial conditions used in 3 -f 1-dimensional cases, where our lattice units are reduced from 
512^ to 512^ In the top panels of Fig. [13] (GRV-M Model in the left panels and GAU-M Model in the 
right panels), we illustrate the evolution of zero- mode and the variances of cr, and in the bottom panels we 
plot the energy density (at t = 3.5 x 10^ in the left-bottom panel and at t = 1.7 x 10^ in the right-bottom 
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panel) instead of the charge density to compare with the Q-ball profiles at zero-temperature, which we 
obtained in FIG. 3 and FIG. 7 in [l^. The colour bars in the bottom panels of Fig. [T3l illustrate the 
values of energy density. Note that we are using the same parameters for GRV-M Model as the ones used 
in [l^ , whilst the potential for GAU-M Model used there is a generalised version of our present potential 
Eq. pi}|) . so the profiles in GAU-M Model should look similar only qualitatively, but not quantitatively. 
From the top panels, we can also see, in particular GRV-M Model, the scaling exponent evolution during 
the driven turbulence stage after the pre-thermalisation ends as confirmed in the top panels of Fig. 1121 
The subsequent evolution, however, is different each other and also unique apart from a characteristic 
free turbulence stage. These features of thermalisation process are caused by stable nonlinear solutions, 
namely "Q-balls" ; in GRV-M Model (left panels) , the variance does not evolve that much after the driven 
turbulence stage ends and we can see thin walled like charged lumps in the end, see the left-bottom panel. 
In GAU-M Model (right panels) the variance has a step-like evolution, at which stage we confirmed that 
two (or sometimes more) charged lumps collide and merge into a larger lump. The collision rate is very 
low since the motions of these "heavy" bubbles are nonrelativistic, but we expect that there will be only 
one single Q-ball left ultimately as similar as the GRV-M case. Generally, we observe that almost all of 
the total energy is trapped into these lumps, where we also confirm that the total charge is absorbed 
into these lumps, as reported in [s^, [3l|. As the "thin- wall" Q-balls in GAU-M Model do not have an 
extremely thin- wall thickness p^ . the profiles seen in the right bottom panel do not have such a thin- 
shell thickness. Note that the "thick-wall" Q-balls in GAU-M Model may suffer from classical instability 
and fission against spatial perturbations around the Q-ball solutions, and decay into smaller Q-balls as 
opposed to the case of "thick-wall" Q-balls in GRV-M Model. The reader should also notice that the 
potential for GAU-M Model in the present case is different from Eqs. (40) and (41) in [l^, which may 
change the classical stability of the Q-balls in the "thick-wall" limit. Furthermore, the stability of Q-balls 
is related to their own charge Q so that the initial ratio, E/{mQ), can also cause the different evolution. 
Therefore, we believe that the evolution is very sensitive to the parameters of models used and the initial 
conditions. It is worth mentioning, in the left-bottom panel, that the value of charge density within 
the charged cluster is slightly larger than the value of the thin- wall Q-balls in the zero-temperature case 
(compare to right bottom panel of FIG. 3 in ^1^). We believe that this is because this charged cluster 
appears in the thermal background, in which the thermal effects change their profiles. 

Let us recap our findings in this section. We have shown in both GRV-M and GAU-M models that 
the AD condensate that has a negative pressure is generally unstable against linear fiuctuations, and the 
fluctuations evolve exponentially. The condition for the presence of the negative pressure corresponds to 
the existence condition of Q-balls, and under our initial conditions shown in TABLE [H we observed that 
almost all of the total charge is trapped into a single (and a few) spherical lump(s) ("thermal Q-balls") 
in the end of our numerical simulations. In the intermidiate regions between the initial exponential 
amplification stage and thermalisation stage in the presence of the nonlinear solutions, we identified that 
the driven turbulence is active; we then found the scaling exponent evolution for the variance of a, and 
we saw that this stage lasts much longer than the case of tachyonic reheating. 

IV. CONCLUSION AND DISCUSSION 

In this paper, we discussed both analytically and numerically two main issues: Affleck-Dine (AD) dy- 
namics and their subsequent nonequilibrium dynamics in the presence of nonlinear solutions. We showed 
that the AD dynamics has the same features as orbital motions of planets, replacing the gravitation force 
by an isotropic harmonic oscillator force. As the relativistic correction on the Newtonian potential gives a 
precession for the planetary orbit motion, the orbits of AD fields are disturbed by the nonrenormalisable 
and quantum correction terms. Note that the essential origin of these corrections is physically different. 
In the presence of a negative pressure of the AD condensate, we have shown that the condensate is 
classically unstable, and the evolution of the system is similar to the dynamics of reheating the Universe, 
i.e. pre-thermalisation, bubble collisions and thermalisation. We found that the thermalisation process 
occurs in the presence of charged lumps, which merge into a single (or a few) "thin-walled (5-ball(s)" 
absorbing most of the homogeneous charge distributed initially in a lattice. 

In Sec. ini we introduced two phenomenological models which are motivated by the MSSM, i.e. gravity- 
mediated (GRV-M) model and gauge-mediated (GAU-M) model. We obtained the frequencies of the 
rotation for the nearly circular orbits, and showed that the condensate can have a negative pressure in 
both cases, see Sec. Ill Al Furthermore, we checked numerically our analytic results with the various cases 
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FIG. 13: (Color online) Left panels (GRV-M Model) and right panels (GAU-M Model) in 2 + 1 dimensions: the 
top panels show the evolution of zero-mode (red-solid lines) and the variations for o (dotted-dashed purple hues). 
In the bottom panels, we plot the energy density (at t = 3.5 x 10^ in the left-bottom panel and at t = 1.7 x 10^ 
in the right-bottom panel) instead of the charge densty to compare the Q-ball profiles seen in FIG. 3 and FIG. 7 
of [l^ where the colour bars illustrate the values of energy density. We can see that almost all of the charge is 
trapped into bubbles which may be "thin-wall" Q-balls, recall that we are imposing a periodic boundary condition. 



in both a non-expanding and expanding universe. 

Our analytic expressions have a number of advantages. In the existing literature on preheating for 
complex scalar fields [13, , the motion of the complex scalar field is assumed to be of an elliptical form, 
but their ansatz does not hold (compare our expressions in Eqs. (|B5l [BT6| and Eq. ^ and their ansatz). 
In the multi-fiat direction cases, our analytic expressions of the AD field give the exact Mathieu equation 
when the interaction term between the AD field (p and another field Xi which parametrises another flat 
direction, is given by where is a coupling constant between them. The previous literature 

[stI suggested that the resonant SUSY preheating for nearly circular orbits is not effective since the 
characteristic dimensionless quantity q is much less than unity, recalling that broad resonant preheating 
(nonadiabatic evolution) occurs for q ^ 1. This statement also holds for our case when the orbit of the 
AD field is nearly circular since q (x where is the third eccentricity of the orbits, recalling that 
nearly circular orbits correspond to the case of ^ 1. 

We obtained the successful ansatze, Eq. for nearly circular orbits in an expanding universe (see 
also the top panels in Fig. [4]), but our analytical expressions could be improved by the action variable 
technique as a real scalar field case [s^. These issues on understanding the analytic orbit forms are 
related to the dynamics of spinning scalar fields, which can be responsible for the early- and late- time 
exponential expansions of the Universe (spinfiation [bO] and spintessence 'gi']) since the AD condensate 
can possess a negative pressure, which can satisfy the condition of slow-roll inflation, w < —1/3. It has 
been discussed in [gl] about an oscillating field responsible for dark energy (see a recent review [11]), and 
it gives a constraint on the power of a power- law potential in order to obtain the attractor solutions (6^ . 
As in real scalar fields, a complex scalar field has been investigated, see for example Following our 
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analytical work, one can investigate the further analysis on dark energy for a complex scalar field and 
their late evolution in order to place constraints on parameters of the models, avoiding Q-ball formation. 

In Sec. mil we explored the late evolution of AD fields in Minkowski spacetime in both GRV-M and 
GAU-M models. As a usual nonequilibrium dynamics, we proposed that the dynamics of the Q-ball for- 
mation goes through three distinct regimes: pre-thermalisation, bubble collision (driven turbulence) and 
thermalisation. We showed analytically that the AD condensate is unstable against spatial perturbations 
when the condensate has a negative pressure, and the perturbations grow exponentially. The presence 
of the negative pressure satisfies the existence of Q-balls as well as the fact that the sound wave of the 
perturbation has an imaginary value of the sound speed. Assuming the adiabatic linear evolution, we 
have analytically presented that the perturbations for the most amplified mode fc = fcm in Eq. ((25|) grow 
with the exponent Sm in Eq. (|26p . which we obtained by taking the average over one rotation of the 



orbits of the AD field. In the previous literature [32, l5l|i these values were obtained by ignoring the 
nonrenormalisable term and by assuming that the orbit is circular. By including the nonrenormalisable 
term and considering more general elliptic orbits, we recovered their results as the leading order of our 
solutions in Sec. IIII A 21 We also showed that the nonlinear time is delayed compared to the time which 
the authors [2^ obtained, since the other modes are not well developed when the most amplified mode 
starts to grow exponentially. With our 3 + 1-dimensional numerical lattice simulations that are run for 
much longer time with much larger simulation sizes than the past lattice simulations in [sol. [3ll [3^ [ssl. [Slj . 
our analytic results are well checked. In addition, we found that the adiabatic condition is violated at the 
beginning stage of the linear perturbations as seen in broad resonant preheating. In the driven turbu- 
lence stage, we observed that many bubbles form and collide/merge into larger bubbles in both GRV-M 
and GAU-M models. Note that these bubbles are nothing to do with bubbles due to first order phase 
transition. By concerning the variance of the radial field cr, we have seen that the evolution follows a 
scaling exponent law as a signature of the driven turbulence fs^l ■ As opposed to the case of tachyonic 
preheating, this driven turbulence stage, in our case, lasts for longer time, which may be caused by the 
presence of classical nonlinear solutions, i.e. "Q-balls". We saw in 2 -I- 1-dimensional numerical results 
that a thermalisation stage actually exists where the evolution for the variance of a field has a different 
scaling law from the one which appears in the driven (first) turbulence stage. We believe that quantum 
effects should be non-negligible in this late turbulence stage, and the classical thermalisation process, 
in our case, should be different from the corresponding quantum-mechanical thermalisation. Since the 
thermalisation process is generally extremely long, a lattice simulation in an expanding background en- 
counters serious problems in the ultra-violet limits; thus, we ignored Hubble expansion in this paper. By 
considering the quantum-mechanical effects as well as Hubble expansion, it is worth to investigate the 
cosmological consequences as our future work. 

In the context of a (p)reheating scenario, it has been suggested [13 that the collision of bubbles during 
the driven turbulence stage can be an effective source of gravitational waves, which can be detected by 
LIGO and LISA in the near future. We noticed that this analysis should be applicable to the 
same driven turbulence stage of the Q-ball formation, which was initially proposed in |18<] . The problem 
of gravitational waves emitted in the fragmentation stage has been discussed while the analysis of 
the gravitational wave emissions in the driven turbulence stage still remains. 
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APPENDIX A: PERTURBATIONS ON MULTIPLE SCALAR FIELDS 



In this appendix, we obtain Euler-Lagrange equations for multiple scalar fields t^" with a symmetric 
field space metric Gab(0) — Gbai'p), following the notations (gsI. [69|. Our aim is to obtain equations of 
motion for the background homogeneous (zero- mode) fields (p'^{t) and the perturbed fields 6Lp°'{t,x.) in a 
fixed unperturbed background (Friedmann-Robertson- Walker) metric, g^i/ = diag(— 1, a{t), a{t), a{t)) 
where a{t) is scale factor of the Universe and H — a/a is Hubble parameter. Here, an over-dot denotes 
the time-derivative. As the simplest nontrivial example of the multiple scalar fields, we find equations of 
motion for a complex scalar field (j) = tre*^ where a and 6 are real scalar fields and the system possesses 
a U(l) symmetry. 

Let us start off with the action 



S 



(fx. 



-g 



(Al) 



where g = det (gf^u) and V{ip) is a potential for ip. Because of the action principle, we obtain the 
Euler-Lagrange equation for (p 



1 



--d„ 



-ggp^CbduVl = 25''"^-^-^^^"^-'^ + 



(A2) 



and the energy momentum tensor 

T^^i, = GabOf^ip^'d^ip'' + g^^ 



-yGabdprda^' - V{ip) 



(A3) 



Here, we defined Gab.c 



dG 



and so on. The energy density and pressure can be given by T^j^ [6S 



PE = ^-g^'Gabd^.v'd^v" + V{v), 
p = -ig^''G,b9^^"9,/-y(^). 



(A4) 
(A5) 



By pertubing the fields as ip'^ — (p'^{t) + (5(^°(i,x) where \(p\ 3> \S(p\, the homogeneous part gives, from 
Eq. dH]), 

(A6) 



^95° + SHip" + G°Vfa = 0, 
at 

where the covariant derivative, D/dt, can be defined by the "Christoffel symbols" 7^^ = ^G'^'^x 

{Gdc,b + Gdb,c — Gbc.d)l thus, |j<y3° = ^tp" + IbcV^v'^- ^n the other hand, we can obtain the equations 
of motion for the pertubed fields 5ip from Eq. (|A2p 

D 



(T^ S^- - 7-^^<^V' V + (Vn^d^v'' = G'^'Gbc.dG'^VJip'', (A7) 



where we used "^Sip"- — Stp'^ + j^^ip''6ip'^ , defined the "Riemann tensors" as 7^^^ = 7^^ c ~ 76c d + Icelbd 
7de7bc' ^'^'^ denoted the covariant derivative as a usual notion Notice that we used Vb = ^^(<^) ■ 



Sip" 



_rV_ 



a + 3H& - cr6|2 



When the system has a 0{2) ~ U{1) symmetry as tpa ~ yr, 9j and a fiat field metric is Gab 
diag(l, tJ^), we can obtain 722 — — ct; 7i2 — 721 — ^jo. Then, Eq. (|A6p with a potential Vip) gives 

da ' 

^ ■ 0. (A9) 



(A8) 



e + me + -te 

a 



Here, the third term in Eq. (|A8|) corresponds to the "centrifugal force" due to the spin in a field space, 
and the third term in Eq. (|A9p corresponds to the "CoUiori force" . In addition, the energy density and 
pressure can be given by from Eqs. (jA4[ IA5P 



PE 



0-2 + cr2^2 



a^ + cr2^2 



V. 



(AlO) 
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Furthermore, Eq. (|A7p gives 



da + 3Hda-{[ — ] + - ^ ] da - 2a96e = 0, (All) 
2ct\ • /V\^ „„ 261 



'SH + —j 60 - i—j S9+ — [aSa - aSaj = 0. (A12) 

We use Eqs. (IA8[ IA9[) to concern with the orbits of AD condensates in Sec. |lTl and use Eqs. (jAlU IA12P 
to investigate the hnear spatial instabihty of the condensates in Sec. 



APPENDIX B: THE ORBIT OF AN AFFLECK-DINE "PLANET 



In this appendix, we obtain an exact orbit form in a quadratic potential case when Hubble expansion 
is assumed to be small and adiabatic. The orbit of a AD field, or more precisely an eccentricity of the 
elliptic motion in the complex field-space, is determined by the initial charge and energy density. In 
order to obtain analytic expressions of the orbit in more general potential cases in which we are more 
interested, we restrict ourself to work in Minkowski spacetime and on the orbit which should be nearly 
circular. We then obtain the perturbed orbit equation and necessary conditions for closed orbits, where 
the orbits come back to their original positions after some rotations around the minimum of the effective 
potential. By including the effects of Hubble expansion, we shall introduce ansatze in Sec. IIIB| which 
are inspired by our solutions obtained in Minkowski spacetime, and our numerical results support the 
ansatze, assuming that the rotation frequency W is always much greater than Hubble expansion H [47| . 



1. The exact orbit in an expanding universe 



The exact orbit expressions of an AD field in an expanding universe can be obtained with a quadratic 
potential, 




where M is a mass of the field (j) and we rescaled the field a, a{t) = f j o-{t). From now on, we solve 

the orbit equations, Eq. ([6]), for a{t) at first, and then solve them for u{9), replacing the time-dependence 
by a phase variable 9. We then show that the orbits are of usual elliptic forms with a third eccentricity 
for a{t) and u{9). 



a. The orbit for a{t) 

In this subsection we obtain an expression for the orbit (j{t) by solving Eq. Substituting Eq. (|Bip 
into Eq. ([6]) and ignoring the terms involving and d/a, we obtain 

^_^+m2^ = ^ ^ = 0, (B2) 
o"' at 

where pE = ^ (w)'^ + ^A'Pa'^ + ^ ^ a^^pE, which is approximately conserved. Since ^^(o'^) = 
+ aa = 2/5e — 2M^a^, Eq. (jB2[) leads to a harmonic oscillator form, 

whose solution is 

^^{t) = ^+^cos[2Af(i + to)], (B4) 
= ^(l + e'cos[2M(i + to)]). (B5) 
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Here, B is some constant value and we set to as a time when the AD field starts to rotate. We also defined 
a third eccentricity = = , where the apocentral and pericentral points are, respectively, 

given by ct^^^ = + A and CT^^n = — A. Notice that the circular orbit case corresponds to = 0, 
which implies that ir^^^ = cr'^im and also note that the eccentricity is real and has a value between 
and 1 in the present quadratic potential [7l|. 
We can obtain a period r of this orbit, 

T = — . (B6) 

Substituting Eq. (|B4p into pE, we obtain A = ^ "^j^^^ — —. From the above expressions for and A, we 
can obtain ^-'^'^ — Vl — e^- Using this and Eq. (IB4p . it ends up with 

Pe 



^ a2 " l + s^cos[2M{t + to)Y ^ ' 

For the circular orbits — 0, 9 is time-independent as we can expect, and the ratio, pe/{Mpq), is unity. 
While for the radial orbits £^ = 1, ^ = as we can also expect, and we can find pe/{Mpq) 3> 1. 



b. The orbit for u{e) =d'^{e) 

What follows is that we express a{t) as a function of 9 by using the second expression in Eqs. (pi IB4[) . 
We then obtain 

tan(6l-6'o) = ?^tan(M(i + <o)), (B8) 

^max 

where Oq is an integration constant and we used the following integral formula, J ai+a^ 



a2 cos X 



_2 Arctan (°i-°^)t^"(f ) ] with some real values ai and 02. Without loss of generality, we can 

choose ^0 = ^^0 = 0, which imply that the orbit at i = 0, r/2 gives, respectively, = 0, 7r/2, recalling 
Eq. (1B6|. By comparing Eq. jlll to Eq. (|B8]), we obtain 

~ 2 ^2 

(9) = — , TT^, (B9j 



T , , 1 cos^ 9 sin^ , 



max mm 



2 + fT^ 

maa; ' ^ min / 1 ^2 



max mtn 



e cos 



(20)) . (Bll) 



Hence, we can see that = when a = (Jmax and 9 — 71/2 when a = dmin- Finally, we obtained the 
expressions for the orbits as usual elliptic forms in Eqs. (jB51 IBllj) . For the circular orbits, = 0, we 
can obtain v?' = const, from Eq. (jBlip as we can expect. 



2. The nearly circular orbits in Minkowski spacetime 



Without Hubble expansion, we can investigate a nearly circular bounded orbit of an AD field in 
general potentials which satisfy Eq. Because of the above reasons, we concentrate on non-expanding 
background in this subsection, and obtain a time-dependent expression for the nearly circular orbits as 
Eq. (jBSp . We then find the expression that depends on the phase 9 as Eq. (jBlOp . Moreover we obtain 
conditions for closed orbits, in which the perturbations are expanded up to 1st order (for the complete 
proof of the condition up to 4th order, see Bertrand's theorem ^33]). 
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a. The orbit for a (t) 

In Minkowski spacetime, we can find an expression for the orbit a{t) in a general potential V{a) as 
in Eq. (|B5p . Notice that the tilde variables are the same as un-tilde ones in the present non-expanding 
background. Recall the equation of motion Eq. ^ in Minkowski spacetime, 

a+^ = 0. (B12) 
da 

Suppose that the orbit that is nearly circular as a{t) = acr + S{t), where acr 3> \S\, recalling acr is defined 
by Eq. ([4]). Substituting this expression of cr into Eq. (jB12p and keeping S terms up to 1st order, we 
obtain a harmonic oscillator form 

S + (B13) 

where the reader should recall the condition, Eq. ([5]), for the bound orbits, and W is constant since we 
are working in Minkowski spacetime. 
Thus, the solution of Eq. (IB13P is 

S{t) = acrBcos{Wt), (B14) 

where i? is a small positive dimensionless constant, i.e. B ^ 1 due to acr S> \6\, and we set the 

differentiation constant as to ensure that a{0) — amax- We can find that amax = o'cr(l + B), amin =^ 

acr{l - B), and a,nax<Jmin ^ <Jcr (l + 0{B'^)) ,which give B = IZllZZZ ^ = """'"2"""" ' and 2B ~ 
2 _ 2 

j2'°^_^j"'" = J where we used the definition of the third eccentricity. We can check that the condition, 

2B ~ <c 1, is consistent with the fact that the orbit is nearly circular. Since amax = ^min = and pE 
is constant, we can equate B with pE and pq using Eqs. (|B141 [5]): 



P 2{pE-V+M) ^ 2{pE-V+M) ^,e^ 

where a prime denotes the differentiation with respect to a. Finally, we obtain 

a^{t) = a^cr (1 + cos{Wt) + 0{e'^)) , (B16) 
where W is given by Eq. ([5]) (compare with Eq. (|B5|) ). Now, we can define the period r 

27r 

which reproduces the case for W = 2M in Eq. (|B6p . Using Eqs. ([2l|4]), we can also find 

0~—^ 7 7- (B18) 

1 + £2 COS (Wt) ^ ' 

Using Eq. ([3]), let us compute the pressure of this AD condensate whose orbit is described by Eq. (|B16[) . 

^22 

By expanding (cr) around a ~ a^ and using Eq. (iBlel) . we obtain {a) ~ (a^) + cos (Wt) + 

( W'^ — -^r- ) cos^ (W^i) + . • . , where we assumed that the higher order terms in V- are negligible. 
Therefore, 

p ~ Yll^si^ (1 _ 2 cos^ {Wt)) - V{acr) + f 1 - 2£' cos [Wt) + ^e^ cos^ [Wt) 
^ (p) ^ -V[a„) + -^. (B19) 

Here we defined an averaged value over an one rotation in the orbit, Eq. (jB16[) . namely (X) = i /J^ dtX{t) 
where X is some quasi-periodic quantity and r is determined by Eqs. (|B171 15)). The result, Eq. (|B19|) . 
can be easily understood by the fact that the averaged pressure corresponds to the value at a = acr since 
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the orbit oscillate around acr and (icr = 0, c.f. a real scalar field case [47|. Similarly, we can obtain the 
averaged energy density 

where we can see the contribution from the term involving . Hence, the averaged equation of state is 
given by 

A - , ^ . (B2:) 



b. The orbit for u{e) = (T'^{e) 



In order to obtain a 6'-dependent expression of the orbit as Eq. (IBlOp . let us switch the variable a to 
u{9) = l/<7{9). In Minkowski spacetime, where we can again drop the tilde variables there, the orbit 
equation Eq. ^ is 



d0^ 



I dV _ 
— — = J{u) 
07, du 



(B22) 



Let uo, which is independent of 6, be the value of a circular orbit (i.e. uq = 1/acr)- We then consider 
an orbit u{6) that deviates slightly from uq with a fluctuation ri{0), i.e. u — uq + rj, where uq ^ |?7|. 
Since ^ = = ^^jgr-, Eq. (|B22[) imphes that uq = J{uo). By expanding J{u) around u = uq, we obtain 
J{u) ~ Uq + 7] + . . . , where we are evaluating the differentiations at uq. Hence, we can obtain the 

perturbed orbit equation for ri{9) 



d'^rj 
d9^ 



+ f3^V = 0, 



(B23) 



where /3 = 1 — which should be positive for bounded orbits. Note that this condition, f3 > 0, is 

equivalent to the previous condition, Eq. ([5]), since 



= ^^W^ 

Pq 



W + aV" 



V 



(B24) 



where we used the fact V = &t u = dcr from Eq. (|4]) . The solution of Eq. (|B23|) is 

■q = uqC cos{/36 + 9q), 



(B25) 



where C and Oq are constants, and < C ^ 1 due to the fact that uq 3> |?7|. We can then show C ^ B 
by equating the value of C with pg and pE- Substituting u into pE and expanding V{u) around u = uq 
up to second order, we can find 



2(pB-y+(l/uo)) 



d^V{l/u) 



ph 



= B 

2 ' 



(B26) 



where we used 



dV+{u) 
du 



Uq 



_ dV{u) 
du 



Uq 



PqUq = from Eq. (jl]). The relation, C = A, is obtained by 



changing the variable u back to a (compare Eq. (jB26p with Eq. (|B15[) ). 
Let us choose 6*0 = vr in Eq. ()B25|) . then we obtain 



u 

,,2 



Uo (1 - Ccosipe)) , 

ul [l-2Ccos(/36') +C'(C2)] 



(B27) 
(B28) 



Notice that < C ^ 1 which is consistent with the condition for nearly circular orbits, ^ 1 as we have 
seen in appendix IB 2 al and Eq. (|B26|) . We can also find that a^ax — -ffc ^'^^ (39 = and amin = ji^ 
for (39 = n. 
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To show that the orbit u{6) in Eq. (|B28[) has a similar form as Eq. (jBlip . let us compute the fol- 
lowing relations: cr^„^ + (t^„ ~ 2al^ (l + ©(C^)) , al,^^ - a^^^ ~ 4a^^C (1 + 0{C)) and cr,^„^cr^„ ~ 

2 I 2 

(j^j. (l + ©(C^)). Hence, u§ ~ "^'^j*^ Jj"'" and 2C ~ e^, which imply that Eq. (jB28P is of similar orbit 

form as Eq. (jBlip . As we computed going from Eq. (jBlOp to Eq. (jBlip . where for this case we deduce 
Eq. (jBlip from Eq. (|B10|) . we finally obtain 

„ cos2 |6l sin^ |6l 
u^^ „ ' +^^- B29) 

max ^ min 

In the next subsection, we obtain the conditions for closed orbits using Eq. (|B29p [70| . 



c. Conditions for closed orbits and equations of state 
Let us define an angle (f>, which is the phase difference as the orbit goes from 77 = uqC to ry = —uqC, 

(B30) 




where we used Eq. (jB24|) . For closed orbits, the angle must have the value that is tt multiplied by a rational 
number, i.e. $ = tt^ where g, r S Z; therefore, (3 should be the rational number. In order to obtain 

the cr-independent value for $, potentials can be of the forms, ^^^^-{+const.), In (cr/m^)^ {+const.), 
and etc. Here, M and are constant real values, and we should have I < —2, < I for bound orbits, 
whereas we may have —2 < ^ < for bound orbits when < 0, recalling Eq. ([5]). The constant terms 
in the potentials add an extra energy for the orbits, and it does not play a significant role, so that we 
consider the potentials without the constant terms. The former power-law potential case, V = ^ , 
gives 

which implies that the closed orbits exist for / = (—1), 2, 7, .... Using Eqs. (jB19| |B20|) and Eq. ([5]), we 
obtain 

(P)^^^^^^, (,,)^([±^^, ^.Jil±3M!^^ (B32) 

which implies that the bound orbits of the AD condensate has a negative pressure for I < 2. In the 
computation of pE, Eq. (|B20[) . we safely ignored the term. The bound orbits for I = (—1,) 2 are 
closed. For the quadratic potential case / = 2, the averaged pressure is zero, in which the AD condensate 
corresponds to an example of nonrelativistic cold dark matter In addition, using Eqs. (|B321 [B2T|) 
we can find 

(-) - (B33) 
On the other hand, the latter logarithmic potential case, In (cr/m^)^, leads to 

^-T^i^ ,B34) 

which corresponds to the former power-law case with I — 0. Similarly, using Eqs. (|B191 IB20p and Eq. ([5]), 
we obtain 



(p)c.m;(l-21n^ , (p^,)^m4 l + 21n^), ^ (B35) 

which implies that the AD condensate has a negative pressure for a^r > ni^ exp (i). In the computation 
oi Pe, Eq- (|B20p . we safely ignored the e'^ term again. Using Eqs. (|B211 |B35)) . we also find 



4m j 
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In Eq. (jB35|) . we cannot clearly see the correspondence with the case for ^ = 0, but we can find ~ 
for <^ Gcr and {w) ~ — 1 for ^ cr as the case with I — 0. 

Let us comment on the pressure when the AD orbit is exact radial, which corresponds to the zero-charge 



density case as for real fields (47[. In this case, the field a{t) coherently oscillates around the vacuum 
if the potential follows a power-law, i.e. y cx cr' for Z > 1, and {w) is the same equation as Eq. (|B33p . 
but it gives a negative pressure for 1 < ^ < 2. Note that the lower bound of I ensures to be a coherent 
oscillation for the radially oscillating AD fields and real scalar fields. 

In summary, we have obtained analytically the explicit expressions, Eqs. (|B51 IBlD) . for the orbit of 
AD fields in a quadratic potential under an expanding universe, and approximately obtained the elliptic 
orbit expressions, Eqs. (jB161 |B29)) . for nearly circular orbits in Minkowski spacetime in potentials which 
satisfy the condition Eq. ([5]) for bound orbits. 



[1] M. Dine and A. Kusenko, Rev. Mod. Ph ys. 76 (2004 ) 1 [a xXiv:hep-ph/0303065| ; A. Riotto and M. Trodden, 
Ann. Rev. Nucl. Part. Sci. 49 (1999) 35 arXiv:hep-ph/9901362 . 



[2] L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56 (1997) 3258 [arXiv:hep-ph/9704452 

J. H. Trasclien and R. H. Brandenberger, Pliys. Rev. D 42 (1990) 2491; Y. Shtanov, J. H. Traschen and 

R. H. Brandenberger, Pliys. Rev. D 51 (1995) 5438 arXiv:hep-ph/9407247 . 
[3] A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5 (1967) 32 [JETP Lett. 5 (1967 SOPUA,34,392-393.1991 

UFNAA,161,6L64.1991) 24]. 
[4] I. Affleck and M. Dine, Nucl. Phys. B 249 (1985) 361; M. Dine, L. Randall and S. D. Tliomas, Nucl. Phys. 

B 458 (1996) 291 arXiv:hep-ph/9507453 . 
[5] L. Randall and R. Sundrum, Nucl. Phys. B 557 (1999) 79 "arXiv:hep-th/9810155'. 

[6] A. Tort, C. Farina and O. de M. Ritter 1989 Eur. J. Phys. 10 220-223 doi: 10.1088/0143-0807/10/3/013 



[7] C. P. Burgess, P. Martineau, F. Quevedo and R. Rabadan, JHEP 0306 (2003) 037 arXiv:h ep-th/0303T70| 
[8] J. G. Rosa and J. March-Russell, Phys. Rev. D 77 (2008) 126004 arXiv:0711.0658 [hep-th]]. 
[9] K. Enqvist and J. McDonald, Phys. Lett. B 425 (1998) 309 |arXiv:hep-ph/9711514i . 
[10] H. P. Nilles, Phys. Rept. 110 (1984) 1. 



[11] R. AUahverdi, B. A. Campbell and J. R. Ellis, Nucl. Phys. B 579 (2000) 355 [arXiv:hep-ph/0001122 
[12] J. McDonald, Phys. Rev. D 48 (1993) 2573 [Erratum-ibid. D 49 (1994) 3066]. 
[13] S. R. Coleman, Nucl. Phys. B 262 (1985) 263 [Erratum-ibid. B 269 (1986) 744]. 
[14] K. M. Lee, Phys. Rev. D 50 (1994) 5333 arXiv:hep-ph/9404293 . 
[15] S. R. Coleman, Phys. Rev. D 15 (1977) 2929 [Erratum-ibid. D 16 (1977) 1248]. 
[16] R. Friedberg, T. D. Lee and A. Sirlin, Phys. Rev. D 13 (1976) 2739. 

[17] T. D. Lee a nd Y. Pang, Phys. Rept. 221 (1992) 251; K. Enqvist and A. Mazumdar, Phys. Rept. 380 (2003) 

99 [arXiv:h ep-ph/0209244 . 

[18] M. I. Tsumagari, E. J. Copeland and P. M. Saffln, Phys. Rev. D 78 (2008) 065021 (arXiv:0805.3233 ' [hep-th]] . 

[19] E. J. Copeland and M. I. Tsumagari, arXiv:0905.0125 [hep-th]. 

[20] M. Laine and M. E. Shaposhnikov, Nucl. Phys. B 532 (1998) 376 arXiv:hep-ph/9804237] . 

[21] R. Banerjee and K. Jedamzik, Phys. Lett. B 484 (2000) 278 arXiv:hep-ph/0005031 . 

[22] K. Enqvist and J. McDonald, Nucl. Phys. B 538 (1999) 321 arXiv:hep-ph/9803380 . 

[23] A. G. Cohen, S. R. Coleman, H. Georgi and A. Manohar, Nucl. Phys. B 272 (1986) 301. 

[24] K. Enqvist, S. Kasuya and A. Mazumdar, Phys. Rev. Lett. 89 (2002) 091301 arXiv:hep-p h/0204270| . 

[25] A. Kusenko and M. E. Shaposhnikov, Phys. Lett. B 418 (1998) 46 arXiv:hep-ph/ 9709492j . 

[26] M. Axenides, S. Komineas, L. Perivolaropoulos and M. Floratos, Phys. Rev. D 61 (2000) 085006 

(arXiv:hep-ph/9910388 ; R. Battye and P. Sutcliffe, Nucl. Phys. B 590 (2000) 329 arXiv:hep-th/0003252 ; 

P. Bowcock, D. Foster and P. Sutcliffe, J. Phys. A 42 (2009) 085403 arXiv:0809.3895 [hep-th]]; C. Elphick, 

Phys. Rev. D 55 (1997) 7749. 

[27] E. Radu and M. S. Volkov, Phys. Rept. 468 (2008) 101 'arXiv:0804.1357' [hep-th]]. 

[28] M. S. Volkov and E. Wohnert, Phys. Rev. D 66 (2002) 085003 arXiv:hep-th/020 5157]; L. Campan elli and 
M. Ruggieri, arXiv:0904.4802 [hep-th]. H. Arodz, J. Karkowski and Z. Swierczynski, arXiv:0907.2801 [hep-th]. 
[29] M. Axenides, E. Floratos, S. Komineas and L. Perivolaropoulos, Phys. Rev. Lett. 86 (2001) 4459 
arXiv:hep-ph/0101193 . 

[30] S. Kasuya and M. Kawasaki, Phys. Rev. D 62 (2000) 02351 2 [arXiv:hep-ph/0002285] ; K. Enqvist, A. Jokinen, 
T. Multamaki and I. Vilja, Phy s. Rev. D 63 (2001) 083501 arXiv:hep-ph/0011134| 7T. Multamaki and I. Vilja, 
Phys. Lett. B 535 (2002) 170 arXiv:hep-ph/0203195 



[31] S. Kasuya and M. Kawasaki, Phys. Rev. D 61 (2000) 041301 [arXiv:hep-ph/9909509 



[32] S. Kasuya and M. Kawasaki, Phys. Rev. D 64 (2001) 123515 arXiv:hep-ph/0106119'. 

[33] J. Berges, AIP Conf. Proc. 739 (2005) 3 arXiv:hep-ph /0409233, ; A. Arrizabalaga, J. Smit and A. Tranberg, 



33 



JHEP 0410 (2004) 017 arXiv:hep-ph/0 409177| ; A. Arrizabalaga, J. Smit and A. Tranberg, Phys. Rev. D 72 
(2005) 025014 arXiv:liep-ph/0503287 . 
[34] R. Micha and I. I. Tkachev, Phys. Rev. D 70 (2004) 043538 arXiv:hep-ph/0403101 . 

[35] G. N. Felder, J. Garcia-Bellido, P. B. Greene, L. Kofman, A. D. Linde and I. Tkachev, Phys. Rev. Lett. 87 

(2001) 011601 arXiv:hep-ph/0012142 ; J. Garcia-Bellido, M. Garcia Perez and A. Gonzalez- Arroyo, Phys. 

Rev. D 67 (2003) 103501 arXiv:hep-ph/0208228 . 
[36] G. N. Felder and L. Kofman, Phys. Rev. D 75 (2007) 043518 arXiv:hep-ph /0606256|; J. F. Dufaux, 

A. Bergman, G. N. Felder, L. Kofman and J. P. Uzan, Phys. Rev. D 76 (2007) 123517 ' arXiv:0707.0875l 

[astro-ph]]. 

[37] J. Garcia-Bellido, D. G. Figueroa and A. Sastre, Phys. Rev. D 77 (2008) 043517 arXiv:0707.0839 [hep-ph]]. 
[38] A. Kusenko and A. Mazumdar, Phy s. Rev. Lett. 101 (2008) 211301 ajX iv: 0_8_01.4554 [astro-p h] ] ; A. Kusenko, 

A. Mazumdar and T. Muhamaki, larXiv:0902.2T97] [astro-ph.CO]; J. F. Dufaux, larXiv:09 02.2574 [astro- 

ph.CO]. 

[39] J. Bertrand, C. R. 77, 849 (1873). 

[40] Fradkin D M1965 Am. J. Phys. 33 207-11; 1967 Prog. Theor. Phys. 37 798-812. 

[41] Laplace P S de 1799 Traite de Mecanique Celeste; Runge C 1919 Vektoranalysis VI pg. 70 (Leipzig: S. Hirzel, 
Leipzig); W Pauh 1926 Z. Phyisik 36 336-63; Lenz W 1926 Z. Phyisik 24 197-207; Heintz W H 1974 Am. J. 
Phys. 42.1078-82 See also Goldstein H 1975 Am. J. Phys. 43 737-8; ibid. 1976 44 1123-4. 

[42] P. Stehle, M. Y. Han Phys. Rev. 159 (1967) 1076. 

[43] P. W. Higgs, J. Phys. A 12 (1979) 309. 

[44] H. I. Leemon, J. Phys. A 12, 489 (1979). 

[45] K. Enqvist, A. Jokinen and J. McDonald, Phys. Lett. B 483 (2000) 191 [arXiv:hep-ph /0004050 . 

[46] A. Jokinen, arXiv:hep-ph/0204086 

[47] M. S. Turner, Phys. Rev. D 28 (1983) 1243. 

[48] N. Bevis and M. Hindmarsh, LATfield Web page, 2006, 'http: / /www.latfield.org/ 1 



[49] http : / / www.nottingham.ac . uk / physics / research /particles / people/Mitsuo / welco me. html | 
[50] A. Pawl, Nucl. Phys. B 679 (2004) 231 arXiv:hep-ph/0412033 . 

[51] K. Enqvist, S. Kasuya and A. Mazumdar, Phys. Rev. D 66 (2002) 043505 tarXiv: hep-ph/0206272) . 

[52] M. Grana and E. Calzetta, Phys. Rev. D 65 (2002) 063522 arXiv:hep-ph/011024 4]. 

[53] S. Y. Khlebnikov and I. I. T kachev, Phys. Rev. D 56 (1997) 653 .arXiv:hep-ph/9701423] . 

[54] |http://www.vapor. ucar.edu] 

[55] S. Y. Khlebnikov and 1. I. Tkachev, Phys. Rev. Lett. 77 (1996) 219 arXiv:he p-ph /9603378| . 

[56] E. J. Copeland, E. W. Kolb and K. M. Lee, Phys. Rev. D 38 (1988) 3023. 

[57] Z. Chacko, H. Murayama and M. Perelstein, Phys. Rev. D 68 (2003) 063515 arXiv:hep-ph/0211369 ; R. Al- 

lahverdi and A. Mazumdar, JCA P 0708 (2007) 023 arXiv:h ep-ph/0608296; ; R. AUahverdi and A. Mazumdar, 

Phys. Rev. D 78 (2008) 043511 [arXw:0 802.4430 [hep-ph]]^ 

[58] M. Postma and A. M azumdar, JCAP 0401 (2004) 005 tarXiv:hep-ph /0304246| ; A. Pawl, Phys. Lett. B 581 

(2004) 231 |arXiv:hep-ph/0411363 . 

[59] M. Berkooz,^jrH. Chung and T. Volansky, Phys. Rev. D 73 (2006) 063526 arXiv:hep-ph/0507218,^^ 

[60] D. A. Easson, R. Gregory, D. F. Mota, G. Tasinato and 1. Zavala, JCAP 0802 (2008) 010 arXiv:07 09.2666l 

[hep-th] ] . 

[61] L. A. Boyle, R. R. Caldwell and M. Kamionkowski, Phys. Lett. B 545 (2002) 17 'arXi v:astro-ph/0105318] . 



[62] A. R. Liddle and R. J. Scherrer, Phys. Rev. D 59 (1999) 023509 arXiv:astro-ph/9809272j 

[63] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15 (2006) 1753 arXiv:hep-th/ 0603057| . 
[64] E. J. Copeland, A. R. Liddle and D. Wands, Phys. Rev. D 57 (1998) 4686 arXiv:gr-qc/9711068' . 
[65] S. Kasuya, Phys. Lett. B 515 (2001) 121 arXiv:astro-ph/0105408 ; M. Alimohammadi and H. Mohseni 
Sadjadi, Phys. Rev. D 73 (200 6) 083527 arXiv:hep-th/0602268. ; H. Wei, R. G. Cai and D. F. Zeng, Class. 



Quant. Grav. 22 (2005) 3189 |arXiv:hep-th/0501160 
[66] 'http://www.ligo. calt ech.edu/[ 
[67] http://lisa.nasa.gOv/l 

[68] M. Sasaki and E. D. Stewart, Prog. Theor. Phys. 95 (1996) 71 arXiv:ast ro-ph/9507001 
[69] S. Weinberg, "Cosmology," Oxford, UK: Oxford Univ. Pr. (2008) 593 p 
[70] Whittaker E T 1937 A Treatise on the Analytical Dynamics of Particles and Rigid Bodies 4*'' ed. (Cambridge: 

Cambridge University Press). 
[71] In an inverse-squared central force, the first eccentricity can be larger than equal 1, which corresponds to the 

cases where the orbits are parabola or hyperbola. 



